A Passivity-Based Velocity Control Method of Hardware-in-the-Loop Simulation for Space Robotic Operations

: The hardware-in-the-loop (HIL) simulation is an important approach to test space robotic operations, rendering virtual free-ﬂoating dynamics on robotic facilities. However, this approach suffers from velocity divergence due to intrinsic time delay in the control loop. This paper proposes a passivity-based control strategy to handle the simulation divergence. A HIL simulation system with an industrial robot is modelled and its divergence problems are analyzed through numerical simulations. Then, through representing the HIL simulation system in a passivity network of view, the passivity observer (PO) of the dynamic system is established. The PO includes the effect of a real contact damping on energy ﬂow of the passivity network during a contact process. Thus, a passivity controller is deﬁned. Moreover, a real-time estimation method for contact damping is presented. Finally, collision experiments against a virtual wall and real collision experiments are both implemented. The experimental results show that the simulation divergence due to the time delay can be prevented by the proposed control strategy, and that the velocity characteristics with high ﬁdelity are rendered on the HIL simulation system.


Introduction
Space robotic operations play a key role in future space missions, such as the removal of space debris in orbit and maintenance tasks on defected satellites [1,2]. These robotic operations all require intensive simulations on the ground in order to validate robotic mechanisms and control strategies in space. Thus, it is of utmost importance to reproduce zero-gravity conditions on the ground. Within this context, there are three categories of simulation technologies, namely full numerical simulations, full physical simulations, and hybrid simulations [3]. The full numerical simulation adopts a software to emulate a space dynamic process, which is a low-cost approach. However, it is not suitable for rendering multiple degrees-of-freedom (DOFs) manipulations because its accuracy depends on mathematical models with unknown physical parameters [4]. The second category includes air-bearing-based testbeds [5,6], water-based neutral buoyancy [7], and parabolic flights [8]. These approaches are usually expensive and not flexible on the changes in simulation parameters. The hybrid simulation, which is also referred to as hardware-in-theloop (HIL) simulation, combines hardware tests with software simulation. This approach has both merits of high flexibility and accuracy, being an attractive approach for simulating complicated robotic tasks in space [1,3].
Many HIL simulators have been developed for the verification of space robotic operations. Shimoji et al. [9] proposed the first HIL simulator, which consisted of a multi-DOF docking-mimicking mechanism, a 5-DOF translational target, and simulation software. Mitchell et al. [10] in the Johnson Space Center presented a 6-DOF parallel robotic facility (SDTS) to simulate rendezvous and docking (RvD) processes. Takhashi et al. [11] presented distinguish the damping force from the elastic force in a space contact process. Since the work carried out by the damping force is just the dissipated energy, the practical contact process including energy dissipation can be rendered on the HIL simulator. Moreover, the energy increase due to time delay can be prevented.
The contribution of this paper is to propose a passivity-based control strategy for solving velocity divergence of a HIL simulation system due to time delay. This approach provides the virtual damping of a traditional PC with an exact physical meaning by means of the real-time identification of contact damping. A HIL simulation system with an industrial robot is modelled and its existing problems due to time delay are analyzed through numerical simulations. Then, the passivity control method based on damping estimation is proposed. Finally, the experiments are implemented and the results validate the proposed method. The remainder of this paper is organized as follows. The HIL simulation system is introduced in Section 2. The control method is proposed in Section 3. Section 4 gives the experiment and discussion. Finally, conclusions are given in Section 5.

The HIL Simulation System
To simulate the free-floating dynamics, a HIL simulator using a 6-DOF industrial robot is proposed, as observed in Figure 1. The mechanical system consists of an industrial robot, a facility robot, a guide rail, a pair of docking-mimicking mechanisms, a satellite mockup, and a 6-axis force/torque (F/T) sensor. The industrial robot is mounted on the guide rail to extend the translational range of the simulator. The docking imitation mechanisms are designed as a collision rod and a collision frame. To produce a point contact between the frame and the rod, the wedge-shaped inner sides of the collision frame are fabricated. The facility robot with the collision rod is installed on the end-effector of the industrial robot. The 6-axis F/T sensor is between the collision frame and the satellite mockup. The specifications of the HIL simulation system are shown in Table 1.  The control hardware of the HIL simulation system is shown in Figure 2. A main controller was adopted to run the entire control program. The control software was developed using the TwinCAT software platform, a real-time PC-based control system from Beckhoff company. The main controller communicates with six servo drives of the industrial robot and a six-axis F/T sensor in real time via the EtherCAT industrial network protocol. In addition, the main controller communicates with the drive of the rail and the space robot by means of the RS-422 serial data standard and CANopen protocol, respectively.

Main
Controller ADS The control hardware of the HIL simulation system is shown in Figure 2. A main controller was adopted to run the entire control program. The control software was developed using the TwinCAT software platform, a real-time PC-based control system from Beckhoff company. The main controller communicates with six servo drives of the industrial robot and a six-axis F/T sensor in real time via the EtherCAT industrial network protocol. In addition, the main controller communicates with the drive of the rail and the space robot by means of the RS-422 serial data standard and CANopen protocol, respectively.  Figure 1. The HIL simulator.
The control hardware of the HIL simulation system is shown in Figure 2. A main controller was adopted to run the entire control program. The control software was developed using the TwinCAT software platform, a real-time PC-based control system from Beckhoff company. The main controller communicates with six servo drives of the industrial robot and a six-axis F/T sensor in real time via the EtherCAT industrial network protocol. In addition, the main controller communicates with the drive of the rail and the space robot by means of the RS-422 serial data standard and CANopen protocol, respectively. The control approach for realizing a zero-gravity environment is shown in Figure 3. There is a hardware layer and a software layer. The hardware layer consists of the industrial robot, the space robot, and the 6-axis F/T sensor. The software layer includes a kinematics module, a motion control subsystem, a free-floating dynamics module, and a timedelay compensation module. The free-floating dynamics module receives the collision force from the 6-axis F/T sensor, and then yields the desired motion trajectory of the simulator. After that, the HIL simulator can emulate the satellite motion process in Cartesian space through the kinematics module. However, there is a time delay between the measured force and the position sent to the robot. The time delay comes from two aspects, The control approach for realizing a zero-gravity environment is shown in Figure 3. There is a hardware layer and a software layer. The hardware layer consists of the industrial robot, the space robot, and the 6-axis F/T sensor. The software layer includes a kinematics module, a motion control subsystem, a free-floating dynamics module, and a time-delay compensation module. The free-floating dynamics module receives the collision force from the 6-axis F/T sensor, and then yields the desired motion trajectory of the simulator. After that, the HIL simulator can emulate the satellite motion process in Cartesian space through the kinematics module. However, there is a time delay between the measured force and the position sent to the robot. The time delay comes from two aspects, namely a dynamic response delay from the industrial robot and a fixed delay from the measurement system. If there is no compensation control, the measured force is directly sent to the free-floating dynamics module, as observed in the blue flow chart of Figure 3. As aforementioned, the time delay will lead to simulation divergence [3,17]. Therefore, it is necessary to propose a compensation module. The kinematics module and the free-floating dynamics module are illustrated in this section and the detailed compensation module will be given in Section 3. namely a dynamic response delay from the industrial robot and a fixed delay from the measurement system. If there is no compensation control, the measured force is directly sent to the free-floating dynamics module, as observed in the blue flow chart of Figure 3. As aforementioned, the time delay will lead to simulation divergence [3,17]. Therefore, it is necessary to propose a compensation module. The kinematics module and the freefloating dynamics module are illustrated in this section and the detailed compensation module will be given in Section 3.

Kinematics of the Simulator
The coordinate systems of the HIL simulator are established, as observed in Figure  4. {Obase} is the global coordinate system; {OIO} and {OIE} are the base coordinate system and the end-effector coordinate system of the industrial robot, respectively; {OME} is the assembly coordinate system on the simulator of the servicing satellite; {OF} is the assembly coordinate system on the mockup of the target satellite; {OS} is the coordinate system of the F/T sensor.

Kinematics of the Simulator
The coordinate systems of the HIL simulator are established, as observed in Figure 4. {O base } is the global coordinate system; {O IO } and {O IE } are the base coordinate system and the end-effector coordinate system of the industrial robot, respectively; {O ME } is the assembly coordinate system on the simulator of the servicing satellite; {O F } is the assembly coordinate system on the mockup of the target satellite; {O S } is the coordinate system of the F/T sensor. measurement system. If there is no compensation control, the measured force is directly sent to the free-floating dynamics module, as observed in the blue flow chart of Figure 3. As aforementioned, the time delay will lead to simulation divergence [3,17]. Therefore, it is necessary to propose a compensation module. The kinematics module and the freefloating dynamics module are illustrated in this section and the detailed compensation module will be given in Section 3.

Kinematics of the Simulator
The coordinate systems of the HIL simulator are established, as observed in Figure  4. {Obase} is the global coordinate system; {OIO} and {OIE} are the base coordinate system and the end-effector coordinate system of the industrial robot, respectively; {OME} is the assembly coordinate system on the simulator of the servicing satellite; {OF} is the assembly coordinate system on the mockup of the target satellite; {OS} is the coordinate system of the F/T sensor.   where s and c denote the sine and cosine functions. Therefore, the forward kinematics of the robot is described as where the coordinate system {O 0 -x 0 y 0 z 0 } and the {O IO } coincide; {O 6 -x 6 y 6 z 6 } and {O IE } coincide. Since the industrial robot includes three consecutive axes that intersect at a point, it is well-known that this kind of robot has a closed-form solution. The detailed inverse kinematics can be found in [29].

Free-Floating Dynamics
Space contact dynamics is a complex full 3-D collision physical process. Compared with the satellites, docking mechanisms or robotic operation tools have significantly smaller sizes. The physical contact that results in robotic operations can be regarded as a point contact. Thus, there is only contact force between the docking mechanisms. Without considering the friction, the contact force occurs along the normal of the contact surface. Accordingly, the full collision can be regarded as a 1-D collision along the resultant force at the contact point. Figure 5 shows coordinate systems of two free-floating satellites. {O so } is the global coordinate of the satellite system, which corresponds to the {O base } in Figure 1. {O ss } and {O st } are the coordinate systems at the centers of mass (CoMs) of two satellites, respectively. The relative relationship between the {O ss } and the {O st } systems corresponds to the {O IE } in Figure 4. {O fs } is the coordinate sytem of the F/T sensor, which corresponds to {O S }. In addition, force analysis during space robotic operations is also depicted in Figure 3. After the space robot on the servicing satellite operates the target satellite, there exists a pair of contact force, F ps and F pt . The contact force can be measured by the F/T sensor. Thus, the force and moment applied on the center of mass (CoM) of the satellite, F ps and M ps , F pt and M pt , can be obtained. Substituting them into the following free-floating dynamic equations leads to the desired motion trajectories of the satellites.  Figure 5. Coordinate systems and force analysis for free-floating satellites.
The translational dynamics of a free-floating satellite can be denoted by and the rotational dynamics is written as are the forces and moments applied on the center of mass (CoM) of the satellite, respectively; are the linear veloc- The translational dynamics of a free-floating satellite can be denoted by Aerospace 2022, 9, 368 7 of 18 and the rotational dynamics is written as where F C (t) ∈ R 3×1 and τ C (t) ∈ R 3×1 are the forces and moments applied on the center of mass (CoM) of the satellite, respectively; v(t) ∈ R 3×1 and ω(t) ∈ R 3×1 are the linear velocity of the CoM and the angular velocity of the satellite; M is the mass of the satellite and I ∈ R 3×3 is the inertia matrix of the satellite, with respect to the coordinate system at the CoM. The 6-axis T/T sensor is usually adopted to measure the contact force and moment. Thus, to describe the force wrench in regard to the CoM frame, the following transformations are given by where c s R is the rotation matrix between the sensor frame and the CoM frame; r c,s is the vector from the origin of the sensor frame to the origin of the CoM frame; F s and τ s are the measured force and moment after measurement compensation, respectively.

Passive Network
The traditional control system view of a HIL simulation system consists of the trajectory generator, the controller/actuator, the sensor, and the plant, as observed in Figure 6a. Since there is a response delay for the HIL simulator, the robot interacts with the environment with a delayed velocity v 2 . After each collision, the trajectory generator also obtains a delayed force, owing to the measurement delay of the 6-axis F/T sensor. According to the passivity theory, a traditional control view can be analyzed in terms of energy flow through representing it in a network point of view, as depicted in Figure 6b. The virtual energy flow can be described using a conjugate pair, such as force and velocity, voltage and current. Energy is written as the integral of inner product between the conjugate input and output [25]. Note that the 'Energy' may or may not correspond to a real physical energy. The main idea of a passivity-based control is to observe the energy flow through a passivity observer (PO) in real time and to dissipate the energy produced owing to the time delay through a passivity controller (PC).
Aerospace 2022, 9, x FOR PEER REVIEW 8 of 19 energy flow can be described using a conjugate pair, such as force and velocity, voltage and current. Energy is written as the integral of inner product between the conjugate input and output [25]. Note that the 'Energy' may or may not correspond to a real physical energy. The main idea of a passivity-based control is to observe the energy flow through a passivity observer (PO) in real time and to dissipate the energy produced owing to the time delay through a passivity controller (PC).

PO and PC
One must consider a dynamic system ℒ with one power port Figure 7), which energetically interacts with the environment. F(τ) and v(τ) are force-like and velocity-like variables, respectively. Note that the conjugate variables in the HIL control system are discrete-time values.

PO and PC
One must consider a dynamic system L with one power port [F(τ), v(τ)] ∈ R n × R n (right-side of Figure 7), which energetically interacts with the environment. F(τ) and v(τ) are force-like and velocity-like variables, respectively. Note that the conjugate variables [F(τ), v(τ)] in the HIL control system are discrete-time values.
delay F comp F act generator (b)

PO and PC
One must consider a dynamic system ℒ with one power port [ ( ), ( )] Figure 7), which energetically interacts with the environment. F(τ) and v(τ) are force-like and velocity-like variables, respectively. Note that the conjugate variables in the HIL control system are discrete-time values.
One-port with time delay According to definition of the passivity [25], ℒ with initial energy storage E(0) is passive if and only if, where Tsam is the sampling time. According to definition of the passivity [25], L with initial energy storage E(0) is passive if and only if, where T sam is the sampling time. Equation (7) means that the energy applied to a passive system must be positive all the time. If E(k) < 0, L is producing energy and then the dynamic system becomes unstable. Moreover, the amount of the produced energy is −E(k). Since the exact value of energy generated is known, a time-varying element can be defined to dissipate the energy produced. This variable element is called a PC. There are serial or parallel configurations for the PC design. An admittance PC requires the dynamics of the industrial robot. However, the dynamics of an industrial robot is usually unknown. Thus, a series configuration with an impedance causality is adopted for the PC on the left-side of Figure 7.
Because there is a time delay, the robotic simulator interacts with the environment with a delayed velocity v des (k − µ), which produces a force v des (k − µ)Z e . Here, Z e is a damping-like variable. In addition, there is also contact damping dissipation between the imitation docking mechanisms. Accordingly, the energy flow can be calculated and the PO can be defined.
The PO of L is the energy flow through the power port [F(k), v(k)] at each time-step, which is given by where F ek (k) is the elastic contact force; F s (k) is the force after the measurement compensation, defined in Equation (10); F α is the PC compensation force; c d is the contact damping.
The identification method of c d will be given in the following subsection. In addition to the time delay, the discrete-integration degrades the behavior of the free-floating dynamics [28]. The work carried out by a rebound force during the acceleration process is larger than the work carried out by a resistance force during the deceleration process, owing to the mismatching of the integration steps of two stages. Therefore, to decrease the effect of the discrete-time integration on velocity divergence, Equation (9) calculates separately the elastic forces during the acceleration and deceleration stages. As aforementioned, there is a fixed time delay for a 6-axis F/T sensor. The time delay can be described by a first-order model with a delay time τ m as follows: where F mea (k) and M mea (k) are the force and moment measured by the sensor; G(s) = 1 + τ m s and L −1 [G(s)] is the inverse Laplace transformation. For a multi-DoF dynamic system, the energy produced needs be dissipated in the Cartesian space of a robot. Therefore, the energy E obs (k) is decomposed into three coordinate axis directions according to the measured force, which is written as For multi-DOF contact in space, since the directions of v des and F ek are not the same, the projection of the desired contact velocity to the contact force is used to calculate the PC compensation force, which is denoted by where α(k) ∈ R 3×1 is the time-varying damping matrix with a diagonal form, and α(k) = diag[α 11 (k), α 22 (k), α 33 (k)]. Evidently, if E obs,i (k) < 0, to dissipate the energy generated at each step, the components of the damping matrix are given by Similarly, since the moment measured by the sensor can be calculated by the contact force at the contact point, the PC compensation moment is written as where r fsp is the vector from the origin of the sensor coordinate system to the contact point. Accordingly, the force and moment compensated by the PC, i.e., the inputs of the space dynamics, are denoted by

Damping Estimation
The above subsection provides a series PC method with impedance causality for preventing the velocity divergence of a HIL simulation. For the impedance configuration, v des (k) = v act (k) is the input and the PC compensation force, F comp (k) and M comp (k) are the outputs. In an ideal case (i.e., the contact damping in Equation (9) is not considered), contact velocities always remain unchangeable during the whole HIL simulation. However, this result does not represent the practical case in which the damping dissipation exists inevitably. Therefore, to reproduce the practical contact dynamics in space, contact damping is necessary to be identified in the HIL simulation.
The position and velocity of the contact point with respect to the global coordinate system in Figure 8 are written as venting the velocity divergence of a HIL simulation. For the impedance configuration, vdes(k) = vact(k) is the input and the PC compensation force, Fcomp(k) and Mcomp(k) are the outputs. In an ideal case (i.e., the contact damping in Equation (9) is not considered), contact velocities always remain unchangeable during the whole HIL simulation. However, this result does not represent the practical case in which the damping dissipation exists inevitably. Therefore, to reproduce the practical contact dynamics in space, contact damping is necessary to be identified in the HIL simulation.
The position and velocity of the contact point with respect to the global coordinate system in Figure 8 are written as  At the time k, given the force Fs(k), the deformation and velocity at the contact point along the direction of contact force can be obtained as Since the sampling time for the force measurement is very short, contact stiffness and contact damping can be regarded as fixed values in one sampling time. Thus, the following relationship can be obtained: where kd and cd are the equivalent contact stiffness and damping along the contact force.
To effectively evaluate the identification algorithm, the facility robot always remains a fixed configuration during the contact process. Therefore, the identified values can be re- At the time k, given the force F s (k), the deformation and velocity at the contact point along the direction of contact force can be obtained as Since the sampling time for the force measurement is very short, contact stiffness and contact damping can be regarded as fixed values in one sampling time. Thus, the following relationship can be obtained: where k d and c d are the equivalent contact stiffness and damping along the contact force.
To effectively evaluate the identification algorithm, the facility robot always remains a fixed configuration during the contact process. Therefore, the identified values can be regarded as an overall stiffness and damping, including the facility robot and the docking-mimicking mechanism. For the above contact process, the Sage-Husa adaptive Kalman filter (AKF) is used to predict and update the estimation values of contact parameters, which is written as where X(k) = [k d (k), c d (k)] T is the state vector; Z(k) = ∆F s (k) is the measurement vector; W(k − 1) is the measurement noise at k − 1; V(k) is observation noise; A(k) is the statetransition matrix from k − 1 to k; G(k) is the system noise matrix; and H(k) = [∆p s (k), ∆v s (k)] T is the observation matrix. Since there are two state variables, both A and G are 2 × 2 identity matrices. The detailed process of solving AKF equations can be found in [30]. Accordingly, the contact damping can be identified when a groups of increments ∆F s (k), ∆p s (k), and ∆v s (k) are given. Although the contact stiffness is simultaneously estimated, the presented control strategy does not require the contact stiffness. In Section 4, the comparison of the identified stiffness and the theoretical value is carried out only to verify the identification method.

Control Strategy
As shown in Figure 9, the whole process of the proposed passivity-based control strategy for the HIL simulation is given as follows: (1) Given the initial displacements, velocities, and accelerations of two satellites, s 0 , v 0 , and a 0 , respectively, the robotic simulator follows the motion trajectory to realize the first collision between the docking imitation mechanisms. (2) The six-axis F/T sensor measures the contact force and moment, F mea and M mea .
Then, the pure time delay τ m caused by the measurement system is compensated by a low-pass filter, and the actual measuring force and moment F s and M s are obtained. (3) By substituting F s , M s , p O p , and v O p into Equations (16)- (22), the contact damping, c d , is identified using the AKF method. (4) According to the identified contact damping, the elastic contact force, F ek , is calculated.
Then, by substituting F ek , M s , v des , and ω des into the PO yields the time-varying damping matrix, and thus the PC compensation force and moment, F α and M α , are calculated through Equations (12) and (14). Accordingly, the compensated force and moment, F comp and M comp , are obtained using Equation (15).

Control Strategy
As shown in Figure 9, the whole process of the proposed passivity-based control strategy for the HIL simulation is given as follows:  (1) Given the initial displacements, velocities, and accelerations of two satellites, s0, v0, and a0, respectively, the robotic simulator follows the motion trajectory to realize the first collision between the docking imitation mechanisms. is identified using the AKF method. (4) According to the identified contact damping, the elastic contact force, Fek, is calculated. Then, by substituting Fek, Ms, des v , and des  into the PO yields the time-varying damping matrix, and thus the PC compensation force and moment, Fα and Mα, are calculated through Equations (12) and (14). Accordingly, the compensated force and moment, Fcomp and Mcomp, are obtained using Equation (15).

Experiment and Discussion
To verify the proposed control strategy, two groups of experiments are carried out on the HIL simulation system, namely collision experiments against a virtual wall, and the real contact experiments between the collision rod and the frame. On one hand, for the virtual collision, contact forces are generated through colliding against a virtual wall, which is modelled as a spring-damper system. Both the contact stiffness and damping of a virtual wall are the fixed values. Through comparing the estimated stiffness and damping with them, the contact parameter estimation method can be verified. On the other hand, for the real contact experiments, the whole control strategy can be verified. Figure  10 shows an example of one real contact process. Since the velocity remains unchangeable

Experiment and Discussion
To verify the proposed control strategy, two groups of experiments are carried out on the HIL simulation system, namely collision experiments against a virtual wall, and the real contact experiments between the collision rod and the frame. On one hand, for the virtual collision, contact forces are generated through colliding against a virtual wall, which is modelled as a spring-damper system. Both the contact stiffness and damping of a virtual wall are the fixed values. Through comparing the estimated stiffness and damping with them, the contact parameter estimation method can be verified. On the other hand, for the real contact experiments, the whole control strategy can be verified. Figure 10 shows an example of one real contact process. Since the velocity remains unchangeable between the front frame and the rear frame, these data are removed in the following experimental curves and the left data are connected to show only the contact process more clearly. between the front frame and the rear frame, these data are removed in the following experimental curves and the left data are connected to show only the contact process more clearly.

Initial velocity
Docking imitation mechanism (Collision rod and frame) Front frame

Collisions against a Virtual Wall
In The virtual wall is a spring-damper system, which is denoted by  The virtual wall is a spring-damper system, which is denoted by where k g = 70 N/mm is the contact stiffness. Since contact damping is a variable with respect to the depth of penetration, a Step function is defined to describe the damping, given by Step(p, 0, 0, p max , c g ) = where c g = 0.1 Ns/mm; p is the depth of penetration; p max is the maximum penetration depth. Figure 11 presents the estimates of contact stiffness and contact stiffness. It is found that contact stiffness and damping are both identified very well. The average of identified stiffnesses is 70.23 N/mm and the maximum of stiffness errors is less than 7%. The average damping is 0.068 Ns/mm, which is close to the given value. The errors of the identified damping are larger than the stiffness errors. Figure 12 gives contact forces and contact velocities. It can be observed that the practical contact forces (i.e., the measured contact forces) remain stable. The maximum of the contact forces after each collision almost remains unchangeable. As mentioned above, the intrinsic time delay of an industrial robot leads to energy increase and simulation divergence. However, the rebound velocities for each collision almost correspond to the desired velocities, which validates the proposed passivity-based control strategy.

Collisions between Docking Imitation Mechanisms
As an example, the masses of the servicing satellite and the target satellite are given as 450 kg and 4000 kg, respectively. The inertia parameters of two satellites and the initial contact velocity are the same as those in collision experiments against a virtual wall. Figure 13a shows the identified contact stiffness. To describe clearly the contact process, a contact frequency can be defined using an average stiffness, as follows:

Collisions between Docking Imitation Mechanisms
As an example, the masses of the servicing satellite and the target satellite are given as 450 kg and 4000 kg, respectively. The inertia parameters of two satellites and the initial contact velocity are the same as those in collision experiments against a virtual wall. Figure 13a shows the identified contact stiffness. To describe clearly the contact process, a contact frequency can be defined using an average stiffness, as follows:

Collisions between Docking Imitation Mechanisms
As an example, the masses of the servicing satellite and the target satellite are given as 450 kg and 4000 kg, respectively. The inertia parameters of two satellites and the initial contact velocity are the same as those in collision experiments against a virtual wall. Figure 13a shows the identified contact stiffness. To describe clearly the contact process, a contact frequency can be defined using an average stiffness, as follows: where m e = m s m t /(m s + m t ) is the relative mass; m s and m t are the masses of the service and target satellites, respectively; k d is the average contact stiffness. The average stiffness for m s = 450 kg is 78.63 N/mm; thus, the contact frequency is obtained as 2.2 Hz. To further validate the control strategy, different contact frequencies are chosen by changing the mass of the servicing satellite. Figure 13b gives the average and the standard deviation of the stiffness estimates with respect to different contact frequencies.
Here, the standard deviation is denoted by where k di is the identified contact stiffness every time. It is found that the averages of the identified stiffnesses for three contact frequencies are very close. In addition, Figure 14 shows the identified contact damping at f = 2.2 Hz and the standard deviation of the damping estimates. It is observed that the averages of the damping estimates for different frequencies are also very close. Note that there exists obvious negative damping during the experiments. The main reason is that force noises lead to large deflections of the estimation error covariance of the AKF. In the control loop, substituting the negative damping factor into the force compensation formula yields the increase in the system energy. Accordingly, the estimated results for different contact frequencies verify the proposed method for contact parameter identification. shows the identified contact damping at f = 2.2 Hz and the standard deviation of the damping estimates. It is observed that the averages of the damping estimates for different frequencies are also very close. Note that there exists obvious negative damping during the experiments. The main reason is that force noises lead to large deflections of the estimation error covariance of the AKF. In the control loop, substituting the negative damping factor into the force compensation formula yields the increase in the system energy. Accordingly, the estimated results for different contact frequencies verify the proposed method for contact parameter identification.  Figure 15 shows the contact velocities and the contact forces at f = 2.2 Hz. There is significant velocity divergence without control because of the time delay. However, the simulation divergence can be prevented by the proposed control strategy. Although there is a slight overshoot of contact velocity for the first collision, the contact velocity converges gradually after two collisions. The time delay in the HIL simulation system leads to a significant force divergence without control, as observed in Figure 15b. To avoid the damage of the imitation docking mechanisms, the experiments without control have to be stopped at the second collision. However, using the proposed passivity-based control strategy, large impact forces at the second collision are prohibited. The contact forces decrease with the decrease in the contact velocities. The contact force convergence occurs. The energy change at f = 2.2 Hz in Figure 16 also shows that the energies without control increase  Figure 15 shows the contact velocities and the contact forces at f = 2.2 Hz. There is significant velocity divergence without control because of the time delay. However, the simulation divergence can be prevented by the proposed control strategy. Although there is a slight overshoot of contact velocity for the first collision, the contact velocity converges gradually after two collisions. The time delay in the HIL simulation system leads to a significant force divergence without control, as observed in Figure 15b. To avoid the damage of the imitation docking mechanisms, the experiments without control have to be stopped at the second collision. However, using the proposed passivity-based control strategy, large impact forces at the second collision are prohibited. The contact forces decrease with the decrease in the contact velocities. The contact force convergence occurs.
The energy change at f = 2.2 Hz in Figure 16 also shows that the energies without control increase significantly. Using the proposed control method, the energy decreases little by little due to the damping dissipation.  To further illustrate the effect of the control strategy on the simulation divergence, a coefficient of velocity (CoV) and a coefficient of contact force (CoF) are defined. The CoV is defined as the ratio of the rebound velocity to the desired velocity (i.e., the initial velocity) [17]. During practical contact experiments, both the theoretical contact stiffness and damping are unknown. Thus, the desired contact force cannot be obtained in practical experiments. Therefore, the CoF is defined as the ratio of the maximum contact force during each collision to the maximum contact force during the first collision. Figure 17 shows the CoVs and CoFs during the experiments with respect to different contact frequencies.
It is found that the convergence phenomena occur more obviously with the increase in contact frequency. Meanwhile, contact forces for different contact frequencies show similar changes to the contact velocities. Accordingly, the proposed control strategy can deal with the simulation divergence for different contact frequencies. To further illustrate the effect of the control strategy on the simulation divergence, a coefficient of velocity (CoV) and a coefficient of contact force (CoF) are defined. The CoV is defined as the ratio of the rebound velocity to the desired velocity (i.e., the initial velocity) [17]. During practical contact experiments, both the theoretical contact stiffness and damping are unknown. Thus, the desired contact force cannot be obtained in practical experiments. Therefore, the CoF is defined as the ratio of the maximum contact force during each collision to the maximum contact force during the first collision. Figure 17 shows the CoVs and CoFs during the experiments with respect to different contact frequencies.
It is found that the convergence phenomena occur more obviously with the increase in contact frequency. Meanwhile, contact forces for different contact frequencies show similar changes to the contact velocities. Accordingly, the proposed control strategy can deal with the simulation divergence for different contact frequencies.

Conclusions
This study proposes a passivity-based velocity control method for HIL simulation systems with time delay. A HIL simulation system is presented using a common industrial robot and its existing problems are investigated. It is found that there is obvious velocity divergence without delay compensation, even if a very small time delay exists. After that, the control strategy is conducted for the velocity divergence, which consists of the PO, the PC, and the damping identification module. The energy flow of the HIL simulation system is analyzed and then the PO is defined. The contact damping identified in real time is introduced into the PC. Thus, the modified damping factor is capable of representing the energy dissipation due to the damping force during a real physical contact process in space. Finally, collision experiments against a virtual wall and collision experiments between docking-mimicking mechanisms are both implemented. The experimental results show that the velocity divergence due to time delay has been prevented by the proposed control strategy, and that the velocity characteristics with high fidelity are obtained.
Currently, only the velocity characteristics of the HIL simulation system are involved in the control strategy. However, the contact force characteristics cannot be handled by the proposed control method, which is determined by the passivity principle. The passivity-based control adjusts the contact velocities to conform to the energy conservation; thus, the accuracy of the contact force is inevitably sacrificed. Moreover, the pure velocity control is likely to produce large contact forces, which will damage the tested instruments. Therefore, a hybrid force/velocity control method for the HIL simulation requires further study.

Conclusions
This study proposes a passivity-based velocity control method for HIL simulation systems with time delay. A HIL simulation system is presented using a common industrial robot and its existing problems are investigated. It is found that there is obvious velocity divergence without delay compensation, even if a very small time delay exists. After that, the control strategy is conducted for the velocity divergence, which consists of the PO, the PC, and the damping identification module. The energy flow of the HIL simulation system is analyzed and then the PO is defined. The contact damping identified in real time is introduced into the PC. Thus, the modified damping factor is capable of representing the energy dissipation due to the damping force during a real physical contact process in space. Finally, collision experiments against a virtual wall and collision experiments between docking-mimicking mechanisms are both implemented. The experimental results show that the velocity divergence due to time delay has been prevented by the proposed control strategy, and that the velocity characteristics with high fidelity are obtained.
Currently, only the velocity characteristics of the HIL simulation system are involved in the control strategy. However, the contact force characteristics cannot be handled by the proposed control method, which is determined by the passivity principle. The passivitybased control adjusts the contact velocities to conform to the energy conservation; thus, the accuracy of the contact force is inevitably sacrificed. Moreover, the pure velocity control is likely to produce large contact forces, which will damage the tested instruments. Therefore, a hybrid force/velocity control method for the HIL simulation requires further study. Data Availability Statement: Data will be made available upon reasonable request by the corresponding author.

Conflicts of Interest:
The authors declare no conflict of interest.