1. Introduction

2. Effect of Different Delays on SITS Stability

Figure 1. Ground semi-physical test process of SITS.

Figure 1. Ground semi-physical test process of SITS.

Figure 2. The 1-DoF undamped free vibration of SIT.

Figure 2. The 1-DoF undamped free vibration of SIT.

3. Mathematical Modeling and Stability Analysis of Hybrid Simulator

3.1. Mathematical Modeling of Hybrid Simulator

To solve the issues of simulation distortion and system divergence, a hybrid simulator is constructed by changing the semi-physical test architecture. The motion signals in the real system are processed and compensated to the dynamic model of one-dimensional undamped vibration of non-cooperative targets, thereby improving the stability of the system. The actual position signal of SITS is processed and compensated into the dynamic model of SIT to improve the stability of the system.

The dynamic equations of the ground semi-physical experimental system, which takes into account FMD, ECC, and SRD before compensation, are as follows

$$\left\{\begin{array}{l}{F^{\prime }}_x\left(t\right)=m_{\mathrm{e}}\ddot{X}\left(t\right)\\ m_{\mathrm{e}}=\frac{m_{\mathrm{c}}{m}_{\mathrm{s}}}{m_{\mathrm{c}}+m_{\mathrm{s}}}\\ X\left(t\right)\left[1\left(t\right)-1\left(t-T_{\mathrm{c}}\right)\right]=X_{\mathrm{c}}\left(t\right)\\ a_{\mathrm{L}{n}_{\mathrm{L}}}{X}_{\mathrm{R}}^{\left(n_{\mathrm{L}}\right)}\left(t\right)+\cdots +a_{\mathrm{L}1}{X^{\prime }}_{\mathrm{R}}\left(t\right)+a_{\mathrm{L}0}{X}_{\mathrm{R}}\left(t\right)=b_{{\mathrm{Lm}}_{\mathrm{L}}}{X}_{\mathrm{c}}^{\left(m_{\mathrm{L}}\right)}\left(t\right)+\cdots +b_{\mathrm{L}1}{X^{\prime }}_{\mathrm{c}}\left(t\right)+b_{\mathrm{L}0}{X}_{\mathrm{c}}\left(t\right)\\ F_x\left(t\right)=-K_{\mathrm{tc}}{X}_{\mathrm{R}}\left(t\right)-B_{\mathrm{tc}}{X}_{\mathrm{R}}\left(t\right)\\ {F^{\prime }}_x\left(t\right)=F_x\left(t-{\tau }_{\mathrm{s}}\right)\end{array}\right.$$

(19)

where ${F^{\prime }}_x\left(t\right)$ is the input force of the dynamic model, $m_{\mathrm{e}}$ is equivalent mass, $m_{\mathrm{c}}$ is mass of tracking spacecraft, $X_{\mathrm{c}}\left(t\right)$ is displacement command input to SITS considering control cycle, $X_{\mathrm{R}}\left(t\right)$ is actual displacement output acquired from SITS controller, $n_{\mathrm{L}}$ and $m_{\mathrm{L}}$ are respectively the orders of SITS dynamic model, and $n_{\mathrm{L}}\ge m_{\mathrm{L}}$. $a_{\mathrm{L}{n}_{\mathrm{L}}}$ and $b_{{\mathrm{Lm}}_{\mathrm{L}}}$ are respectively the differential equation coefficients, $B_{\mathrm{tc}}$ is the contact damping coefficient.

Overlay the force measurement delay signal with the disposed velocity signal to correct the input of the dynamic model. Let the compensation force be $F_{\mathrm{h}}\left(t\right)$, then rewrite the 6th formula in Equation (19) as follows

$$\left\{\begin{array}{l}{F^{\prime }}_x\left(t\right)=F_x\left(t-{\tau }_{\mathrm{s}}\right)+F_{\mathrm{h}}\left(t\right)\\ F_{\mathrm{h}}\left(t\right)=K_{\mathrm{h}}\left(t\right){\dot{X}}_{\mathrm{R}}\left(t\right)\end{array}\right.$$

(20)

where $K_{\mathrm{h}}\left(t\right)$ is the gain, which can be a fixed value or a time-varying value.

Perform Laplace transform on Equation (19) and Equation (20), and the corresponding control block diagram of the hybrid simulator is shown in Figure 3.

Figure 3. Control block diagram of the hybrid simulator.

Figure 3. Control block diagram of the hybrid simulator.

The expressions of $G_{\mathrm{d}}\left(s\right)$ and $G_{\mathrm{f}}\left(s\right)$ in Figure 3 are

$$\left\{\begin{array}{c}{G}_{\mathrm{d}}\left(s\right)=\frac{1}{m_e{s}^2}\\ G_{\mathrm{f}}\left(s\right)=-\left(K_{\mathrm{tc}}+B_{\mathrm{tc}}s\right)\end{array}\right.$$

(21)

Due to the fact that $e^{-\tau s}$ is an infinite dimensional function, Pade approximation is used for linearization, which has good performance even when delay time is large and Taylor series does not converge [27,28,29]. The first-order Pade approximation expressions for $G_c\left(s\right)$ and $G_{\mathrm{s}}\left(s\right)$ are as follows

$$\left\{\begin{array}{l}{G}_c\left(s\right)=\frac{2T_{\mathrm{c}}}{2+T_{\mathrm{c}}s}\\ G_{\mathrm{s}}\left(s\right)=\frac{2-{\tau }_{\mathrm{s}}s}{2+{\tau }_{\mathrm{s}}s}\end{array}\right.$$

(22)

3.2. Model Identification of SITS

$G_R\left(s\right)$ is transfer function of 1D motion of SITS, which mainly includes the forward kinematics model, control model of each joint, and inverse kinematics model, as shown in Figure 4. In practice, $G_R\left(s\right)$ is generally obtained through system identification.

Figure 4. Main constituent models of $G_R\left(s\right)$

Figure 4. Main constituent models of $G_R\left(s\right)$

In order to comprehensively evaluate the performance of simulator, a certain load is configured, and it is expected that as many axes as possible move simultaneously during the identification. In addition, the motion of simulator should be obvious, making it convenient for real-time measurement by devices such as the laser tracker. Set the 1D motion trajectory as P0→P1→P0→P2→P0→P1→..., as shown in Figure 5.

Figure 5. The 1D motion trajectory for model identification of SITS.

Figure 5. The 1D motion trajectory for model identification of SITS.

The sweep path equation of the simulator is

$${\overline{X}}_{\mathrm{c}}\left(t\right)={\overline{A}}_{\mathrm{R}}\sin \left({\overline{\omega }}_{\mathrm{R}}t\right)$$

(23)

where ${\overline{A}}_{\mathrm{R}}$ and ${\overline{\omega }}_{\mathrm{R}}$ are amplitude and angular frequency, respectively.

Control point $P_{\mathrm{v}}$ of simulator is located at P0 at initial moment, and the initial velocity is

$$\left\{\begin{array}{l}{V}_{\mathrm{R}0}={\overline{A}}_{\mathrm{R}}{\overline{\omega }}_{\mathrm{R}}{e}_{\mathrm{P}}\\ e_{\mathrm{P}}=\frac{\overrightarrow{P_0{P}_1}}{\left|\overrightarrow{P_0{P}_1}\right|}\end{array}\right.$$

(24)

Let the coordinates of $P_{\mathrm{v}}$ in simulator base coordinate system be $P_{\mathrm{v}}={\left[\begin{array}{ccc}{p}_x& p_y& p_z\end{array}\right]}^{\mathrm{T}}$, yields

$${\overline{X}}_{\mathrm{c}}\left(t\right)e_{\mathrm{P}}=\overrightarrow{P_0{P}_{\mathrm{v}}}$$

(25)

The joint angular displacement $q_i\left(i=1, 2,\cdots , 6\right)$ corresponding to above Cartesian trajectory can be obtained through robot kinematics model, as follwings

$$\left\{\begin{array}{l}{}_6{}^{0}T{}_s{}^{6}T={\left[\begin{array}{cc}{P}_{\mathrm{v}}^{\mathrm{T}}& 1\end{array}\right]}^{\mathrm{T}}\\ {}_s{}^{6}T=\left[\begin{array}{cc}{I}_{3\times 3}& {}^6\mathbf{P}_0\\ 0& 1\end{array}\right]\\ {}_6{}^{0}T=f\left(q_1,q_2,q_3,q_4,q_5,q_6,a_1,a_2,a_3,d_1,d_4,{\alpha }_1,{\alpha }_3,{\alpha }_4,{\alpha }_5\right)\end{array}\right.$$

(26)

where ${}^6\mathbf{P}_0$ is the coordinate vector of p0 in the sixth axis coordinate system of SITS, ${}_6{}^{0}T$ is homogeneous transformation matrix of SITS.

The Z-transform function $G_R\left(z\right)$ of SITS can be acquired by applying least squares method to identify the above trajectory in the discrete time domain.

The ARX model [30,31,32] of $G_R\left(z\right)$ can be written as

$$\left\{\begin{array}{c}A\left(z^{-1}\right)X_{\mathrm{R}}\left(k\right)=B\left(z^{-1}\right)X_{\mathrm{c}}\left(k-d\right)+v\left(k\right)\\ A\left(z^{-1}\right)=1+a_1{z}^{-1}+a_2{z}^{-2}+\cdots +a_n{z}^{-n_{\mathrm{R}}}\\ B\left(z^{-1}\right)=b_0+b_1{z}^{-1}+b_2{z}^{-2}+\cdots +b_m{z}^{-m_{\mathrm{R}}}\end{array}\right.$$

(27)

where $v\left(k\right)$ is systematical error.

The difference equation corresponding to Equation (27) is

$$\left\{\begin{array}{l}{X}_{\mathrm{R}}\left(k\right)={\phi }_{\mathrm{R}}^{\mathrm{T}}\left(k\right){\theta }_{\mathrm{R}}+v\left(k\right)\\ {\phi }_{\mathrm{R}}\left(k\right)=\left[-X_{\mathrm{R}}\left(k-1\right),-X_{\mathrm{R}}\left(k-2\right),\cdots ,-a_n{X}_{\mathrm{R}}\left(k-n_R\right),\right.\\ {\left.X_{\mathrm{c}}\left(k\right),X_{\mathrm{c}}\left(k-1\right),X_{\mathrm{c}}\left(k-2\right)，\cdots ,X_{\mathrm{c}}\left(k-m_R\right),v\left(k\right)\right]}^{\mathrm{T}}\\ {\theta }_{\mathrm{R}}={\left[a_1,a_2,\cdots ,a_n,b_0,b_1,\cdots ,b_m\right]}^{\mathrm{T}}\end{array}\right.$$

(28)

If the number of observed samples is $N_{\mathrm{R}}$, the least squares estimate ${\widehat{\theta }}_{\mathrm{R}}$ of ${\theta }_{\mathrm{R}}$ is

$${\widehat{\theta }}_{\mathrm{R}}={\left({\Phi }_{\mathrm{R}}^{\mathrm{T}}{\Phi }_{\mathrm{R}}\right)}^{-1}{\Phi }_{\mathrm{R}}^{\mathrm{T}}{X}_{\mathrm{R}}$$

(29)

where

$$\left\{\begin{array}{l}{X}_{\mathrm{R}}={\left[X_{\mathrm{R}}\left(2\right),X_{\mathrm{R}}\left(3\right),\cdots ,X_{\mathrm{R}}\left(N_{\mathrm{R}}\right)\right]}^{\mathrm{T}}\\ {\Phi }_{\mathrm{R}}={\left[{\phi }_{\mathrm{R}}^{\mathrm{T}}\left(2\right),{\phi }_{\mathrm{R}}^{\mathrm{T}}\left(3\right),\cdots ,{\phi }_{\mathrm{R}}^{\mathrm{T}}\left(N_{\mathrm{R}}\right)\right]}^{\mathrm{T}}\end{array}\right.$$

(30)

3.3. Stability and Steady-State Error Analysis

The control block diagram in Figure 3 satisfies

$$\left\{\begin{array}{c}{F}_{\mathrm{in}}\left(s\right)G_{\mathrm{ds}}\left(s\right)+{F^{\prime }}_x\left(s\right)G_{\mathrm{ds}}\left(s\right)+\frac{{F^{\prime }}_x\left(s\right)}{G_{\mathrm{f}}\left(s\right)G_{\mathrm{s}}\left(s\right)}{F}_{\mathrm{h}}\left(s\right)G_{\mathrm{ds}}\left(s\right)={F^{\prime }}_x\left(s\right)\\ G_{\mathrm{ds}}\left(s\right)=G_{\mathrm{d}}\left(s\right)G_c\left(s\right)G_R\left(s\right)G_{\mathrm{f}}\left(s\right)G_{\mathrm{s}}\left(s\right)\end{array}\right.$$

(31)

The zero and pole form of identified $G_R\left(s\right)$ is given by

$$G_R\left(s\right)=K_R\frac{\prod _{j=1}^{m_{\mathrm{L}}}\left(s-{z_r}_i\right)}{\prod _{i=1}^{n_{\mathrm{L}}}\left(s-{p_r}_i\right)}$$

(32)

where $K_R$ is gain, $m_{\mathrm{L}}$ and $n_{\mathrm{L}}$ are respectively the number of zeros ${z_r}_i$ and poles ${p_r}_i$.

Transfer function $G_{\mathrm{w}}\left(s\right)$ of SITS can be obtained by performing Laplace transform on the second formula in Equation (20) and substituting it into Equation (31) along with Equation (21), Equation (22) and Equation (32).

$$\begin{array}{rl}& G_w(s)=\frac{F_x^{’}(s)}{F_{\text{in}}(s)}\\ & =\frac{\frac{1}{m_e{s}^2}\frac{2T_c}{2+T_{c}s}{K}_R\frac{\underset{j=1}{\overset{m_L}{\Pi }}(s-z_{ri})}{\underset{i=1}{\overset{r_L}{\Pi }}(s-p_{ri})}\frac{2-{\tau }_{s}s}{2+{\tau }_{s}s}[-(K_{tc}+B_{\mathrm{t}\mathrm{c}}s)]}{1-\frac{1}{m_e{s}^2}\frac{2T_c}{2+T_{c}s}{K}_R\frac{\underset{j=1}{\overset{m_L}{\Pi }}(s-z_{ri})}{\underset{i=1}{\Pi }(s-p_{ri})}\frac{2-{\tau }_{s}s}{2+{\tau }_{s}s}[-(K_{\mathrm{t}\mathrm{c}}+B_{\mathrm{t}\mathrm{c}}s)]-K_{\mathrm{h}}s\frac{1}{m_e{s}^2}\frac{2T_c}{2+T_{c}s}{K}_R\frac{\underset{i=1}{\overset{m_L}{\Pi }}(s-z_{ri})}{\underset{i=1}{\overset{r_L}{\Pi }}(s-p_{ri})}}\\ & \end{array}$$

(33)

The characteristic equation of Equation (33) is given by

$$\begin{array}{l}{m}_e{s}^2\left(2+T_{\mathrm{c}}s\right)\left(\prod _{i=1}^{n_{\mathrm{L}}}\left(s-{p_r}_i\right)\right)\left(2+{\tau }_{\mathrm{s}}s\right)+2T_{\mathrm{c}}{K}_R\prod _{j=1}^{m_{\mathrm{L}}}\left(s-{z_r}_i\right)\left(2-{\tau }_{\mathrm{s}}s\right)\left(K_{\mathrm{tc}}+B_{\mathrm{tc}}s\right)\\ -2T_{\mathrm{c}}{K}_{\mathrm{h}}{K}_{R}s\left(2+{\tau }_{\mathrm{s}}s\right)\prod _{j=1}^{m_{\mathrm{L}}}\left(s-{z_r}_i\right)=0\end{array}$$

(34)

From identification results in Section 5, there are two poles and one zero of $G_R\left(s\right)$, so $n_{\mathrm{L}}=2$ and $m_{\mathrm{L}}=1$. Equation (34) is further simplified as

$$\begin{array}{l}{m}_e{s}^2\left(2+T_{\mathrm{c}}s\right)\left(s-{p_r}_1\right)\left(s-{p_r}_2\right)\left(2+{\tau }_{\mathrm{s}}s\right)\\ +2T_{\mathrm{c}}{K}_R\left(s-{z_r}_1\right)\left(2-{\tau }_{\mathrm{s}}s\right)\left(K_{\mathrm{tc}}+B_{\mathrm{tc}}s\right)-2T_{\mathrm{c}}{K}_{\mathrm{h}}{K}_{R}s\left(2+{\tau }_{\mathrm{s}}s\right)\left(s-{z_r}_1\right)=0\end{array}$$

(35)

Define the coefficients of system characteristic equations are $a_{\mathrm{r}0},a_{\mathrm{r}1},\cdots ,a_{\mathrm{r}6}$, then Equation (35) can be rewritten as

$$a_{\mathrm{r}6}{s}^6+a_{\mathrm{r}5}{s}^5+a_{\mathrm{r}4}{s}^4+a_{\mathrm{r}3}{s}^3+a_{\mathrm{r}2}{s}^2+a_{\mathrm{r}1}s+a_{\mathrm{r}0}=0$$

(36)

where

$$\{\begin{array}{l}{a}_{\mathrm{r}0}=-4T_c{K}_R{K}_{\mathrm{t}\mathrm{c}}{z}_{r1}\\ a_{\mathrm{r}1}=2T_{\mathrm{c}}{K}_R[(2+{\tau }_{\mathrm{s}}{z}_{r1})K_{\mathrm{t}\mathrm{c}}-2z_{r1}{B}_{\mathrm{t}\mathrm{c}}]+4T_{\mathrm{c}}{K}_{\mathrm{h}}{K}_R{z}_{\mathrm{r}1}\\ a_{\mathrm{r}2}=4p_{r1}{p}_{r2}{m}_e+2T_c{K}_R[(2+{\tau }_{\mathrm{s}}{z}_{r1})B_{\mathrm{t}}-{\tau }_{\mathrm{s}}{K}_{\mathrm{t}}]-2T_{\mathrm{c}}{K}_{\mathrm{h}}{K}_R(2-{\tau }_{\mathrm{s}}{z}_{r1})\\ a_{\mathrm{r}3}=[2({\tau }_{\mathrm{s}}+T_c)p_{r1}{p}_{r2}-4(p_{r1}+p_{r2})]m_e-2{\tau }_{\mathrm{s}}{T}_c{K}_R{B}_{\mathrm{t}\mathrm{c}}-2T_c{K}_{\mathrm{h}}{K}_R{\tau }_{\mathrm{s}}\\ a_{\mathrm{r}4}=[T_c{\tau }_s{p}_{r1}{p}_{r2}-2({\tau }_s+T_c)(p_{r1}+p_{r2})]m_c\\ a_{\mathrm{r}5}=[2({\tau }_s+T_{\mathrm{c}})-T_c{\tau }_s(p_{r1}+p_{r2})]m_c\\ a_{\mathrm{r}6}=m_c{T}_c{\tau }_s\end{array}$$

(37)

Owing to good stability and anti-interference ability of industrial robot, ${z_r}_i<0$ and $\mathrm{Re}\left({p_r}_i\right)<0$. Considering that ${\tau }_{\mathrm{s}}$ and $T_{\mathrm{c}}$ are relatively small, it can be clearly concluded from Equation (37) that $a_{\mathrm{r}0}>0$, $a_{\mathrm{r}4}>0$, $a_{\mathrm{r}5}>0$ and $a_{\mathrm{r}6}>0$.

In terms of Routh stability criterion, Routh array of Equation (36) is given by

(38)

To make sure the stability of SITS, coefficients of the first column of Routh array must be positive, illustrating as in Equation (39).

$$\left\{\begin{array}{l}{A}_{\mathrm{r}1}=\frac{a_{\mathrm{r}4}{a}_{\mathrm{r}5}-a_{\mathrm{r}3}{a}_{\mathrm{r}6}}{a_{\mathrm{r}5}}>0\\ B_{\mathrm{r}1}=\frac{A_{\mathrm{r}1}{a}_{\mathrm{r}3}-A_{\mathrm{r}2}{a}_{\mathrm{r}5}}{A_{\mathrm{r}1}}>0\\ C_{\mathrm{r}1}=\frac{B_{\mathrm{r}1}{A}_{\mathrm{r}2}-B_{\mathrm{r}2}{A}_{\mathrm{r}1}}{B_{\mathrm{r}1}}>0\\ D_{\mathrm{r}1}=\frac{C_{\mathrm{r}1}{B}_{\mathrm{r}2}-a_{\mathrm{r}0}{B}_{\mathrm{r}1}}{C_{\mathrm{r}1}}>0\end{array}\right.$$

(39)

where

$$\left\{\begin{array}{l}{A}_{\mathrm{r}2}=\frac{a_{\mathrm{r}2}{a}_{\mathrm{r}5}-a_{\mathrm{r}1}{a}_{\mathrm{r}6}}{a_{\mathrm{r}5}}\\ B_{\mathrm{r}2}=\frac{A_{\mathrm{r}1}{a}_{\mathrm{r}1}-a_{\mathrm{r}0}{a}_{\mathrm{r}5}}{A_{\mathrm{r}1}}\end{array}\right.$$

(40)

The steady state error of hybrid simulator can be drawn as

$$\begin{array}{cc}\hfill e_{ss}=& \underset{t\to \infty }{\lim }E\left(t\right)\hfill \\ & =\underset{s\to 0}{\lim }sE\left(s\right)\hfill \\ & =\underset{s\to 0}{\lim }s\left(F_{\mathrm{in}}\left(s\right)+{F^{\prime }}_x\left(s\right)+\frac{{F^{\prime }}_x\left(s\right)F_{\mathrm{h}}\left(s\right)}{G_{\mathrm{f}}\left(s\right)G_{\mathrm{s}}\left(s\right)}\right)\hfill \end{array}$$

(41)

When $F_{\mathrm{in}}\left(s\right)$ is a step function with the expression as $F_{\mathrm{in}}\left(s\right)=1/s$, Equation (41) yields

$$e_{ss}=\underset{s\to 0}{\lim }s\left(1+G_{\mathrm{w}}\left(s\right)+\frac{G_{\mathrm{w}}\left(s\right)F_{\mathrm{h}}\left(s\right)}{G_{\mathrm{f}}\left(s\right)G_{\mathrm{s}}\left(s\right)}\right)F_{\mathrm{in}}\left(s\right)=0$$

(42)

From the above analysis, it can be concluded that when appropriate parameters are selected, the hybrid simulator proposed is convergent and can achieve stable ground experiments.

4. Delay Compensation Method of Hybrid Simulator

Figure 6. Figure 6. Delay compensation principle diagram of hybrid simulator.

Figure 6. Figure 6. Delay compensation principle diagram of hybrid simulator.

5. Simulation Study

Table 1. The parameters and corresponding values required for numerical simulations.

5.1. Simulation Verification of Delays Effect on the SITS Stability of 1D Undamped Vibration

To show the performances of system instability under three different delays FMD, ECC, and SRD, the discrete method is used for solution. The simulation step size for all simulations is set to 0.002s, and the response lag model of simulator adopts the transfer function identified in the next simulation.

The comparisons of velocity $V$, force $F$, and power $P$ of SIT under different delay factors are shown in Figure 7. Figure 7(a) shows the comparison of speed and force between idea case and only considering FMD. Figure 7(b) shows the speed and force comparison between FMD and ECC. Figure 7(c) shows the speed and force contrast between idea case and SRD. Figure 7(d) shows the power comparison of SIT under the conditions of no delay and different delays.

Figure 7. Figure 7. The comparison of velocity, force, and power under different delay factors.

Figure 7. Figure 7. The comparison of velocity, force, and power under different delay factors.

It can be seen from Figure 7 that when there are three different types of delay mentioned above, the 1D undamped system of SIT diverges, but the velocity, force, and power are finite within a specific time interval, which also verifies the analysis and proof of system stability by the aforementioned lags. In the simulation, the delay of FMD has a greater impact than that of ECC, and SRD causes the system to diverge the slowest. Furthermore, the vibration amplitude and energy of SITS under different delay conditions continue to increase compared to the previous cycle. If corresponding compensation is not timely carried out, it will lead to the loss of credibility of ground capture simulation test results, and even damage to the simulator.

5.2. Identification Simulation of SITS

Discretize the sine sweep frequency curves corresponding to the relevant parameters in Table1 and utilize them as control instructions for identifying the discrete domain transfer function of SIST. Then, map the instructions to the Cartesian coordinate system and the joint space of robot through kinematics model. The spatial trajectory of control point Pv is shown in Figure 8.

Figure 8. The spatial trajectory of P_v.

Figure 8. The spatial trajectory of P_v.

Taking the above frequency sweep signal as input and the simulated trajectory as output, $G_R\left(z\right)$ can be obtained through the mentioned identification algorithm, and $G_R\left(s\right)$ can be obtained through transformation. The identified trajectory error $v\left(t\right)$ is shown in Figure 9(a), and its logarithmic frequency characteristic curve is drawn in Figure 9(b).

Figure 9. Figure 9. Model identification results of SITS.

Figure 9. Figure 9. Model identification results of SITS.

As shown in Figure 9, 1D sweep motion of simulator has a large stability margin. Its amplitude frequency characteristics first increase and then decrease with frequency, while the phase lag continues to increase. At 10Hz, its phase lag is −2.87°.

5.3. Simulation Verification of the Effect for SCVG

To verify the effectiveness of the tracking differentiator and optimize its tuning parameters, Figure 10(a) presents the original displacement curve and direct differential signal of the simulator. Figure 10(b) demonstrates differential signal estimation curves of the original signal corresponding to different tracking differentiator parameters. It can be seen that from Figure 10(a) there are significant spikes and oscillations of the differential signal, and directly applying it for feedback control will reduce system stability. From the partial enlarged view in Figure 10(b), it can be seen that when r is set to a larger value, the estimated tracking speed of its differential signal increases significantly, but its amplitude overshoot also increases. When h takes a larger value, the differential signal becomes stable, but the signal will be distorted due to excessive filtering. Therefore, the parameters of tracking differentiator are selected as $\left[r h\right]=\left[1000 0.0034\right]$ after comprehensive consideration, which not only has a fast tracking speed, but also a relatively stable and close to the real differential signal.

Figure 10. Verification of the tracking differentiator effectiveness. ((a)Original displacement curve and differential signal of the simulator, (b)Differential signal estimations under different tracking differentiator parameters.).

Figure 10. Verification of the tracking differentiator effectiveness. ((a)Original displacement curve and differential signal of the simulator, (b)Differential signal estimations under different tracking differentiator parameters.).

To attain better dynamic compensation effects, $t_{\mathsf{\epsilon }}$ plays an critical role to improve the ground experiment accuracy of hybrid simulator in addition to parameter tuning of $K_{\epsilon 0}$, $K_{\epsilon 1}$ and $K_{\epsilon 2}$. Fig.11 shows the comparison of $e_{\mathrm{F}}$ when $t_{\epsilon }$ takes different values in SCVG.

From Figure 11, it can be seen that when $t\le t_{\epsilon }$, the linear compensator takes effect, and the maximum force error after compensation is 6.76N. The nonlinear compensator based on energy observer begins to take effect when $t>t_{\epsilon }$, and the peak deviation between the compensated force and the ideal value when $t_{\mathsf{\epsilon }}=36s$ is slightly larger than $t_{\mathsf{\epsilon }}=26s$, with a duration of 1s and a maximum deviation of 5.39%. Afterwards, the force error is smaller than $t_{\mathsf{\epsilon }}=26s$ and $t_{\mathsf{\epsilon }}=46s$, so it has a better compensation effect when $t_{\mathsf{\epsilon }}=36s$ with $e_{\mathrm{F}}$ reducing from 6.76N to 0.09N, which will be used in the following simulations.

Figure 11. The comparison of $e_{\mathrm{F}}$ when $t_{\epsilon }$ takes different values in SCVG.

Figure 11. The comparison of $e_{\mathrm{F}}$ when $t_{\epsilon }$ takes different values in SCVG.

To prove the advantage of proposed compensator, the velocity error comparison of SIT under different compensation methods, including FGC, LGC and SCVG, is shown in Figure 12.

Figure 12. Comparison of SIT velocity error under different compensator.

Figure 12. Comparison of SIT velocity error under different compensator.

In Figure 12, it can be seen that $e_{\mathrm{V}}$ of LGC is better than that of FGC before the switching control takes effect. The maximum $e_{\mathrm{V}}$ of the former is 0.14mm/s, while the latter is 0.22mm/s. If the nonlinear compensator based on energy observer is not carried out, the velocity error amplitude of LGC exceeds that of FGC at 40.69s. After implementing a nonlinear compensator based on the energy observer, $e_{\mathrm{V}}$ is significantly reduced. At the end of the simulation, the amplitude of $e_{\mathrm{V}}$ is only 0.028mm/s adopting SCVG, while that of LGC and FGC are 0.23mm/s and 0.16mm/s, respectively.

Although $e_{\mathrm{V}}$ is relatively small within the simulation time, due to the continuous attenuation of the vibration velocity, the deviation of its vibration recovery coefficient from the ideal value may become larger. The comparison of the changes of $e_{\mathrm{R}}$ for different compensators is shown in Figure 13.

Figure 13. The comparison of $e_{\mathrm{R}}$ for different compensators during simulation.

Figure 13. The comparison of $e_{\mathrm{R}}$ for different compensators during simulation.

From Figure 13, it can be seen that although $e_{\mathrm{R}}$ of LGC is greater than that of LGC in later stage of simulation, it reduces the energy of the system and thus increases its stability. The $e_{\mathrm{R}}$ of SCVG remains between 0.995 and 1.005, and compared with FGC (1.0203) and LGC (0.9756), it gradually approaches 1 with a value of 0.9995, which also verifies the effectiveness of the proposed compensator.

The movement velocity of SIT is closer to the ideal case without delay due to the adoption of SCVG, which also confirms that the dynamic input force approaches the true value. Figure 14 compares the $e_{\mathrm{F}}$ using different compensators

Figure 14. Comparison of $e_{\mathrm{F}}$ under different compensators.

Figure 14. Comparison of $e_{\mathrm{F}}$ under different compensators.

As can be seen from Figure 14, input force errors of three types of compensators are basically the same in initial phase of the simulation, while the differences enlarge gradually as time increases. The value of $e_{\mathrm{F}}$ applying SCVG decreases the fastest, and peak value $e_{\mathrm{F}}$ at the end of simulation is only 0.29N compared with FGC (0.58N) and LGC (0.77N), which verifies the excellent compensation performance of SCVG in long-term hybrid simulation.

The variation of the observed energy of SIT over time adopting SCVG is revealed in Figure 15.

Figure 15. The variation of ${\epsilon }_{\mathrm{h}}$ adopting SCVG.

Figure 15. The variation of ${\epsilon }_{\mathrm{h}}$ adopting SCVG.

It can be seen from Figure 15 that the value of ${\epsilon }_{\mathrm{h}}$ starts from 0, the amplitude increasing to 0.034w within the first vibration period, and decreases to 8.37×10^-4w in the last simulation cycle, which implies that as the vibration period increases, it eventually tends towards 0, indicating that the system converges. The variation of $K_{\mathrm{h}}$ over time in the proposed SCVG during simulation is shown in Figure 16.

Figure 16. The variation of $K_{\mathrm{h}}$ in SCVG during simulation.

Figure 16. The variation of $K_{\mathrm{h}}$ in SCVG during simulation.

Feedback compensation gain $K_{\mathrm{h}}$ varies linearly with time before the switching controller taking effect. As system energy gradually stabilizes in the later phase of simulation, it actually exhibits nonlinear changes, which means that SITS is corrected in a timely manner after the switch time.

Furthermore, the values of $C_{\mathsf{\epsilon }\mathrm{l}}$ and $C_{\mathsf{\epsilon }\mathrm{h}}$ can also be acquired, which are $-111.80$ and $-110.86$, respectively. The SITS system has two characteristic roots with positive real parts when $K_{\mathrm{h}}=0$, and the system is unstable. It can be concluded that the eigenvalues all have negative real parts by substituting $C_{\mathsf{\epsilon }\mathrm{l}}$ and $C_{\mathsf{\epsilon }\mathrm{h}}$ into characteristic Equation (36) and it simultaneously satisfies the stability constraints of Equation (39), indicating that the SITS is asymptotically stable. Therefore, the proposed hybrid simulation system with SCVG is stable.

6. Conclusions

References

1. Keles O F. Telecommunications and Space Debris: Adaptive Regulation Beyond Earth. Telecommun Policy 2023, 47, 102517. [CrossRef]
2. Takeichi, N.; Tachibana, N. A tethered plate satellite as a sweeper of small space debris. Acta Astronaut 2021, 189, 429–436. [CrossRef]
3. Witze, A. The quest to conquer earth's space junk problem. Nature 2018, 561, 24–26. [CrossRef]
4. A. Nanjangud, P.C. Blacker, S. Bandyopadhyay, Y. Gao. Robotics and AI-Enabled On-Orbit Operations With Future Generation of Small Satellites. P Ieee 2018, 106, 429-439. [CrossRef]
5. N. Inaba, M. Oda, M. Asano. Rescuing a Stranded Satellite in Space-Experimental Robotic Capture of Non-Cooperative Satellites. T Jpn Soc Aeronaut S 2006, 48, 213-220. [CrossRef]
6. E. Papadopoulos, F. Aghili, O. Ma, R. Lampariello. Robotic Manipulation and Capture in Space: A Survey. Front Robot Ai 2021, 8, 686723. [CrossRef]
7. T. Rybus. Obstacle avoidance in space robotics: Review of major challenges and proposed solutions. Prog Aerosp Sci 2018, 101, 31-48. [CrossRef]
8. H. G. Hatch, J. E. Pennington, and J. B. Cobb. Dynamic simulation of lunar module docking with Apollo command module in lunar orbit. NASA Tech. Note TN D-3972, 1967, pp. 1–26.
9. Uyama N, Fujii Y, Nagaoka K, et al. Experimental evaluation of contact-impact dynamics between a space robot with a compliant wrist and a free-flying object. International symposium on artificial intelligence, robotics and automation in space, Turin, Italy. 2012.
10. Wei, Y., Yang, X., Xu, Z., Bai, X. Novel ground microgravity experiment system for a spacecraft-manipulator system based on suspension and air-bearing. Aerosp Sci Technol 2023, 141, 108587. [CrossRef]
11. O. Ma and J. Wang. Model order reduction for impact-contact dynamics simulations of flexible manipulators. Robotica 2007, 25, 397–407. [CrossRef]
12. A. Flores-Abad, O. Ma, K. Pham, S. Ulrich. A review of space robotics technologies for on-orbit servicing. Prog Aerosp Sci 2014, 68, 1–26. [CrossRef]
13. A. Yaskevich. Real time math simulation of contact interaction during spacecraft docking and berthing. J. Mech. Eng. Autom 2014, 4, 1–15.
14. Roe F D, Howard R T, Murphy L. Automated rendezvous and capture system development and simulation for NASA. Modeling, Simulation, and Calibration of Space-based Systems. SPIE 2004, 5420: 118-125.
15. Mitchell, J. D., Cryan, S. P., Baker, K., Martin, T., Goode, R., Key, K. W., Chien, C. H. Integrated docking simulation and testing with the Johnson space center six-degree-of-freedom dynamic test system. In Proc. Space Technol. Appl. Int. Forum, 2008, pp. 709–716.
16. Boge T , Ma O .Using advanced industrial robotics for spacecraft Rendezvous and Docking simulation. In 2011 IEEE international conference on robotics and automation, pp. 1-4.
17. O. Ma, A. Flores-Abad, and T. Boge. Use of industrial robots for hardware-in-the-loop simulation of satellite rendezvous and docking. Acta Astronaut 2012, 81, 335–347. [CrossRef]
18. Wei YQ, Yang X, Bai XL, Xu ZG. Hardware-in-the-loop based ground test system for space berthing and docking mechanism of small spacecraft. P I Mech Eng G-J Aer 2023, 237, 3486-3495. [CrossRef]
19. Zebenay, M., Lampariello, R., Boge, T., Choukroun, D. A new contact dynamics model tool for hardware-in-the-loop docking simulation. 2012.
20. Zebenay M, Boge T, Krenn R, et al. Analytical and experimental stability investigation of a hardware-in-the-loop satellite docking simulator. P I Mech Eng G-J Aer 2015, 229, 666-681. [CrossRef]
21. Qi, C., Ren, A., Gao, F., Zhao, X., Wang, Q., Sun, Q. Compensation of velocity divergence caused by dynamic response for hardware-in-the-loop docking simulator. Ieee-Asme T Mech 2016, 22, 422-432. [CrossRef]
22. Stefano M D, Balachandran R, Secchi C. A Passivity-Based Approach for Simulating Satellite Dynamics With Robots: Discrete-Time Integration and Time-Delay Compensation. Ieee T Robot 2020, 36, 189-203. [CrossRef]
23. Rao, L. S., Rao, P. S., Chhabra, I. M., Kumar, L. S. Modeling and Compensation of Time-Delay Effects in HILS of Aerospace Systems. Iete Tech Rev 2022, 39, 375-388. [CrossRef]
24. Osaki K, Konno A, Uchiyama M. Delay time compensation for a hybrid simulator. Adv Robotics 2010, 24, 1081-1098. [CrossRef]
25. Abiko, S., Satake, Y., Jiang, X., Tsujita, T., Uchiyama, M. Delay time compensation based on coefficient of restitution for collision hybrid motion simulator. Adv Robotics 2014, 28, 1177-1188. [CrossRef]
26. Khusainov, D. Y., Diblík, J., Růžičková, M., Lukáčová, J. Representation of a solution of the Cauchy problem for an oscillating system with pure delay. Nonlinear Oscil 2008, 11, 276-285. [CrossRef]
27. Ricardo L G G, Lagomasino G L. Strong asymptotics of multi-level Hermite-Padé polynomials. J Math Anal Appl 2024, 531, 127801. [CrossRef]
28. Đukić, D., Mutavdžić Đukić, R., Reichel, L., Spalević, M. Optimal averaged Padé-type approximants. Electron T Numer Ana 2023, 59, 145-156. [CrossRef]
29. Motoyama Y, Yoshimi K, Otsuki J. Robust analytic continuation combining the advantages of the sparse modeling approach and the Padé approximation. Phys Rev B 2022, 105, 035139. [CrossRef]
30. Mehmood, K., Chaudhary, N. I., Cheema, K. M., Khan, Z. A., Raja, M. A. Z., Milyani, A. H., Alsulami, A. Design of nonlinear marine predator heuristics for hammerstein autoregressive exogenous system identification with key-term separation. Mathematics 2023, 11, 2512. [CrossRef]
31. Beltran-Perez, C., Serrano, A. A., Solís-Rosas, G., Martínez-Jiménez, A., Orozco-Cruz, R., Espinoza-Vázquez, A., Miralrio, A. A general use QSAR-ARX model to predict the corrosion inhibition efficiency of drugs in terms of quantum mechanical descriptors and experimental comparison for lidocaine. Int J Mol Sci 2022, 23, 5086. [CrossRef]
32. Mehmood, K., Chaudhary, N. I., Khan, Z. A., Cheema, K. M., Raja, M. A. Z., Shu, C. M. Novel knacks of chaotic maps with Archimedes optimization paradigm for nonlinear ARX model identification with key term separation. Chaos Soliton Fract 2023, 175, 114028. [CrossRef]
33. Wen, Q., Wang, M., Li, X., Chang, Y. Learning-based design optimization of second-order tracking differentiator with application to missile guidance law. Aerosp Sci Technol 2023, 137, 108302. [CrossRef]
34. Park J H, Park T S, Kim S H. Asymptotically convergent higher-order switching differentiator. Mathematics 2020, 8, 185. [CrossRef]
35. Feng H, Qian Y. A linear differentiator based on the extended dynamics approach. Ieee T Automat Contr 2022, 67, 6962-6967. [CrossRef]
36. Oliveira T R, Estrada A, Fridman L M. Global and exact HOSM differentiator with dynamic gains for output-feedback sliding mode control. Automatica 2017, 81, 156-163. [CrossRef]