Stability Analysis and Delay Compensation for Space Instable Target Simulator

1. Introduction

With the increase in human space activities, the number of non-cooperative space objects such as defunct satellites and spacecraft debris continues to rise, significantly impacting the utilization of orbital resources and the space operational environment [1,2,3]. Active debris removal (ADR) technology involves the use of tools like robotic arms to repair and clean up space debris, effectively reducing the generation of space debris. This technology proves instrumental in mitigating the threats posed by uncontrolled debris in the space environment and alleviating congestion in orbital space [4,5,6,7].

Direct on-orbit application of ADR technology poses significant risks, generally necessitating comprehensive ground tests to address shortcomings in technical methods and validate the stability. Full physical simulation (FPS), semi-physical simulation (SPS), and digital simulation (DS) are currently the three commonly applied ground experiment approaches. Full-sized simulators with complete physical characteristics are typically constructed in FPS, which incurs high costs. In addition, it is difficult to achieve frictionless movement along the direction of gravity in the simulation of 6-DoF spatial weightlessness motion [8,9,10]. DS, not reliant on physical hardware, conducts verification of debris removal processes only through mathematical modeling [11,12,13]. The mathematical model for actual contact process cannot be accurately established, and it also lacks the capability to optimize and improve the real capture tool, presenting certain limitations. SPS combines the advantages of FPS and DS, incorporating vital physical hardware into the test system to reduce the unnecessary costs of hardware [14,15,16,17,18]. Simultaneously, it utilizes mathematical spacecraft models with different mass and inertia characteristics, making it widely applied in ground experiments for space mechanism.

The ground capture SPS system for space instable target (SIT) is a system with high real-time requirements. SITS is a critical device, and its motion performance is directly related to the accuracy of the simulation. Delays have significant effects on the stability and precision of SPS system, leading to invalid test data or even damage of simulator due to instability.

Researches on the reasons and compensation methods for the delay of simulator in SPS system have gained certain progress. M. Zebenay et al. [19] addressed the response time delay of the industrial robot controller in the European Proximity Operations Simulator (EPOS 2.0 c) and the impact of robot high stiffness on the safety and the authenticity of experiment results. A contact collision model suitable for ground verification is proposed, combining impedance control with mechanical passive compliance, which builds the relationship between robot stiffness and test device stiffness. Based on estimating contact stiffness and damping, M. Zebenay introduced compensation into the dynamic model to achieve 1D hybrid simulation in 2015 [20]. By comparing the maximum contact force and velocity recovery coefficient of 1D movement with an air-floating test system, the feasibility of hybrid simulator was verified. Aiming at response delay of the 6-DoF parallel mechanism in ground SPS platform, Chenkun Qi [21] investigated a real-time online identification method for contact stiffness, which was estimated through differential equations based on force, displacement, and velocity within three contact cycles. Adopting this information as the compensation of dynamic input, system divergence was suppressed with achieving a velocity recovery coefficient of 0.9809. To address the issues of system instability and simulation distortion caused by time delay and discrete integration of dynamics in the On Orbit Servicing Simulator loop, Marco De Stefano [22] proposed a passive controller for adjusting admittance dynamics input force and a speed modulation controller for maintaining energy characteristics. The instability suppression effect in satellite dynamics reproduction was verified on industrial robots. LAV Sanhith Rao [23] analyzed the motion divergence caused by the delay in the three-axis rotational semi-physical test platform for aerospace products. By compensating for the turntable angular velocity based on establishing the system transfer function, they successfully solved the undesired oscillations and inaccurate experiment results. The reason for the increase of system energy was identified through analyzing the phase relationship between force and velocity in complex domain by Kohei Osaki [24]. The force compensation model is established and applied to 1D collision experiments, achieving a recovery coefficient of 0.97 with satisfactory results. However, due to neglecting the velocity phase variation caused by force lag, this mean provides an approximate compensation effect. Satoko Abiko [25] presented a force compensation way based on the velocity recovery coefficient. By observing the velocity change online, they identified the compensation coefficient for input force of dynamic model and applied it to the planar 3-DoF collision tests.

There is a considerable risk to identify stiffness and damping in real-time during actual capture process and the feasibility is weak, especially when the contact stiffness is high, that is, the contact time is also rather short. Moreover, relying solely on adjusting the impedance characteristics of manipulator to indirectly extend the contact time and suppress divergence is not a fundamental approach to solve the divergence issue of simulator. The integral lag in dynamics is just one aspect of control cycle delay. Only settling this problem while neglecting other lag impacts leaves the simulator in a divergent state. Similarly, sampling delay is only a part of system lags. The 3-axis turntable conducts angular velocity control, and its angular velocity is derived from the high-precision measurement sensor. The compensation method is unsuitable for the simulator for the industrial robot only precisely executing displacement commands. It is difficult to predict in advance the contact dynamic characteristics of the captured object. Therefore, the velocity recovery coefficient is also unknown, making it impractical to compensate for the input force of the hybrid simulator based on the velocity recovery coefficient.

The 1D movement of SITS is researched according to the principle of minimum envelope. It has been analyzed and proven that all three types of delay in current SPS can cause the 1D motion divergence of SIT in the absence of damping. On this basis, a hybrid simulation control block diagram is proposed, and its transfer function model is derived. A variable gain differential negative feedback with energy observer is presented to compensate for the input force of dynamic model, and the effectiveness of the proposed method is verified through simulation study.

The rest of the paper is organized as follows: effect of different delays on SITS stability and mathematical modeling and stability analysis of hybrid simulator are described in Section 2 and Section 3, respectively; delay compensation method of hybrid simulator and results obtained by a series of simulations are presented in Section 4 and Section 5, respectively; conclusions are drawn in Section 6.

2. Effect of Different Delays on SITS Stability

The simulator, a 6-DoF industrial robot, is used to conduct spatial motion of SIT in the hardware-in-the-loop based ground capture test system. A six-axis force sensor is installed on the end flange of simulator to measure the contact force during the capture process. After delayed sampling, it is input into the dynamic model of SIT. Motion instructions obtained after dynamic calculation are transmitted to host computer and sent to simulator according to a certain control cycle. Simulator controller decomposes the instructions to motor controllers of each joint, driving each joint motor to generate new contact force, and follow the above process in a cyclic manner until the capture experiment is completed. The ground semi-physical test process of SITS is illustrated in Figure 1.

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

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

During above process, simulator is affected by factors such as force measurement delay, discrete control cycle, and response delay, resulting in reduced motion accuracy, decreased system stability, and even divergence. To illustrate impact of lagging link on simulator stability and find the corresponding compensation method, considering the complexity of the system, one-dimensional motion of SIT is generally analyzed. The undamped free vibration model is adopted here as shown in Figure 2.

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

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

In Figure 2, $V_0$ and $m_s$ are respectively initial velocity and mass of SIT, $K_{\mathrm{tc}}$ is the spring stiffness. SIT is theoretically in a free vibration state, and the system is critically stable. The corresponding dynamic equation is given by

$$\left\{\begin{array}{l}\ddot{X}\left(t\right)+{\omega }_n^{2}X\left(t\right)=0\\ X\left(0\right)=0，\dot{X}\left(0\right)=V_0\end{array}\right.$$

(1)

where natural frequency is ${\omega }_n=\sqrt{K_{\mathrm{tc}}/m_s}$.

2.1. Effect of Force Measurement Delay on System Stability

Let force measurement delay (FMD) be ${\tau }_{\mathrm{s}}$, the vibration differential equation of SIT can be written as

$$\ddot{X}\left(t\right)+{\omega }_n^{2}X\left(t-{\tau }_{\mathrm{s}}\right)=0$$

(2)

Perform Laplace transform on Equation (2), and the corresponding characteristic equation is

$$s^2+{\omega }_n^2{\mathrm{e}}^{-s{\tau }_s}=0$$

(3)

Equation (3) is a transcendental equation, which theoretically has countless characteristic roots. By substituting $s=j\omega$ into the above equation, we can obtain

$$-{\omega }^2+{\omega }_n^2\cos \left(\omega {\tau }_s\right)-j{\omega }_n^2\sin \left(\omega {\tau }_s\right)=0$$

(4)

Separating the real and imaginary parts yields

$$\left\{\begin{array}{c}{\omega }_n^2\cos \left(\omega {\tau }_s\right)={\omega }^2\\ {\omega }_n^2\sin \left(\omega {\tau }_s\right)=0\end{array}\right.$$

(5)

By solving Equation (5), considering the symmetry of the eigenvalues and only considering solutions with $\omega >0$, we can get

$$\left\{\begin{array}{l}\omega ={\omega }_n=\sqrt{\frac{K_{\mathrm{tc}}}{m_s}}\\ {\tau }_{s,n}=\frac{2n\pi }{{\omega }_n},n=0,1,2\cdots \end{array}\right.$$

(6)

Obviously, when ${\tau }_s=0$, $s=\pm j{\omega }_n$ is a pair of conjugate intercepts.

The real component of the derivative of the root locus at all intersection points can be expressed as

$${\left.\mathrm{Re}\left(\frac{ds}{d{\tau }_s}\right)\right|}_{\begin{array}{l}s=j\omega \\ {\tau }_s={\tau }_{s,n}\end{array}}=\frac{2{\omega }^2{\omega }_n^2\cos w{\tau }_{s,n}}{{\left[-2\omega \sin \left(\omega {\tau }_{s,n}\right)-{\omega }_n^2{\tau }_{s,n}\right]}^2+{\left[2\omega \cos \left(\omega {\tau }_{s,n}\right)\right]}^2}>0$$

(7)

A pair of root trajectories enter the right side of the imaginary axis in complex plane when ${\tau }_s=0$, and they will not return to the left side when ${\tau }_s\to \infty$. Therefore, the undamped vibration of SIT is unstable for $\forall {\tau }_s>0$.

To further demonstrate the stability of the time-domain solution in finite time [26], for Equation (2), when $\gamma <\delta$, there is

$$‖X\left(t\right)‖<C，t\in \left[0，T\right]$$

(8)

So, Equation (2) is about the stability of $\left\{\delta ,C,{\tau }_s\right\}$ in finite time interval. Where $\left[0，T\right]$ and $\gamma =\max \left\{‖\phi ‖C,‖\dot{\phi }‖C,‖\ddot{\phi }‖C\right\}$.

In summary, the undamped vibration of SIT with force measurement delay is stable in finite time, but as time approaches infinity, the system gradually diverges.

2.2. Effect of Control Cycle on System Stability

Discrete control cycle keeps the current displacement command for a period of time, which is equivalent to performing a zero order hold on the command signal. The impact model of control cycle $T_{\mathrm{c}}$ can be can be described as follows

$$g_c\left(t\right)=1\left(t\right)-1\left(t-T_{\mathrm{c}}\right)$$

(9)

Performing the Laplace transform on Equation (9) yields

$$G_c\left(s\right)=L\left[1\left(t\right)-1\left(t-T_{\mathrm{c}}\right)\right]=\frac{1}{s}-\frac{e^{-T_{\mathrm{c}}s}}{s}$$

(10)

Undamped vibration model of SIT considering only the effect of control cycle (ECC) is

$$\left\{\begin{array}{l}{F}_x=m_s\ddot{X}\left(t\right)\\ F_x=-K_{\mathrm{tc}}{X}_{\mathrm{c}}\left(t\right)\\ X_{\mathrm{c}}\left(t\right)=X\left(t\right)g_c\left(t\right)\end{array}\right.$$

(11)

where $X_{\mathrm{c}}\left(t\right)$ is the displacement command considering ECC.

By substituting Equation (10) into the Laplace transform of Equation (11) and performing the inverse Laplace transform, it can be obtained that

$$\stackrel{⃛}{X}\left(t\right)+{\omega }_n^{2}X\left(t\right)-{\omega }_n^{2}X\left(t-T_{\mathrm{c}}\right)=0$$

(12)

To prove its stability, we introduce the following inference.

For equation $f_D{\left[\begin{array}{cccc}{f}^3& f^2& f& 1\end{array}\right]}^{\mathrm{T}}=0$, where

$$f_D{\left[\begin{array}{cccc}{f}^3& f^2& f& 1\end{array}\right]}^{\mathrm{T}}=0$$

(13)

If $f$ has multiple non negative roots, then its corresponding time-delay system is unstable, the equation of which is

$$\stackrel{⃛}{x}\left(t\right)+p_2\ddot{x}\left(t\right)+p_1\dot{x}\left(t\right)+p_{0}x\left(t\right)+q_1\dot{x}\left(t-\tau \right)+q_{0}x\left(t-\tau \right)=0$$

(14)

We provide proof for the above inference. The corresponding characteristic equation for the third-order time-delay linear system mentioned in Equation (14) is:

$${\lambda }^3+p_2{\lambda }^2+p_1\lambda +p_0+\left(q_1\lambda +q_0\right)e^{-\lambda \tau }=0$$

(15)

On the basis of the stability criterion of the control system, if there are multiple identical zero real part eigenvalues in Equation (15), then its corresponding system is unstable. Let the pure imaginary root be $j\omega \left(\omega >0\right)$, and by substituting it into Equation (15) and separating the real and imaginary components, we can obtain

$$\left\{\begin{array}{l}-{\omega }^3+p_1\omega =q_0\sin \left(\omega \tau \right)-q_1\omega \cos \left(\omega \tau \right)\\ p_2{\omega }^2-p_0=q_0\cos \left(\omega \tau \right)+q_1\omega \sin \left(\omega \tau \right)\end{array}\right.$$

(16)

Squaring and adding the above equations yields

$${\omega }^6+\left(p_2^2-2p_1\right){\omega }^4+\left(p_1^2-2p_0{p}_2-q_1^2\right){\omega }^2+p_0^2-q_0^2=0$$

(17)

Let $f={\omega }^2$, then Equation (17) can be written as $f_D{\left[\begin{array}{cccc}{f}^3& f^2& f& 1\end{array}\right]}^{\mathrm{T}}=0$. Due to the presence of multiple non negative roots in f, i.e. $\omega$ also has multiple positive real roots, the delay system is unstable.

Evidently, Equation (12) is one type of Equation (17). Therefore, it can be concluded that when considering control cycle from the above inference, the undamped vibration system with 1-DoF is of divergence.

2.3. Effect of Simulator Response Delay (SRD) on System Stability

Due to the reaction time of electromechanical system, there is a certain delay in the motion response of SITS, and the corresponding differential equation can be written as

$$\ddot{X}\left(t-{\phi }_{\mathrm{r}}\left(t,\tau \right)\right)+{\omega }_n^{2}X\left(t-{\phi }_{\mathrm{r}}\left(t,\tau \right)\right)=0$$

(18)

where ${\phi }_{\mathrm{r}}\left(t,\tau \right)$ is time-varying delay function.

In terms of the characteristics of the control system, as the motion frequency increases, the phase delay becomes larger, and the response lag reduces the system adaptability to external interference and phase margin. Therefore, 1-DoF undamped movement of SITS is divergent.

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

The hybrid simulation model for ground capture experiments of SIT is elaborated above in continuous time domain. Due to the fact that the disposals in physical systems is based on sampled values, the compensation force $F_{\mathrm{h}}\left(t\right)$ needs to be discretized. Meanwhile, considering the continuous increase of system energy, an energy observer is introduced to dynamically adjust the gain, and the principle diagram of delay compensation is illustrated in Figure 6.

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

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

In Figure 6,${\widehat{\dot{X}}}_{\mathrm{R}}\left(k\right)$ is the estimation of ${\dot{X}}_{\mathrm{R}}\left(k\right)$ obtained by tracking differentiator, ${\epsilon }_{\mathrm{h}}\left(k\right)$ is the observed energy, $K_{\mathrm{h}}\left(k\right)$ is the gain. The k-th force compensation value $F_{\mathrm{h}}\left(k\right)$ can be expressed as

$$\left\{\begin{array}{l}{F}_{\mathrm{h}}\left(k\right)=K_{\mathrm{h}}\left(k\right){\widehat{\dot{X}}}_{\mathrm{R}}\left(k\right)\\ k=\frac{t}{T}\end{array}\right.$$

(43)

where T is system sampling time.

To quantify the performance of proposed compensation method, evaluation indicators of dynamic input force error $e_{\mathrm{F}}$, SIT velocity error $e_{\mathrm{V}}$, and recovery coefficient $e_{\mathrm{R}}$ are introduced, the expressions of which are given by

$$\left\{\begin{array}{l}{e}_{\mathrm{F}}=E\left(t\right)-{F^{\prime }}_x\left(t\right)\\ e_{\mathrm{V}}={\dot{X}}_{\mathrm{com}}\left(t\right)-\dot{X}\left(t\right)\\ e_{\mathrm{R}}=\frac{{\dot{X}}_{\mathrm{com}}^{+}\left(t\right)}{{\dot{X}}^{+}\left(t\right)}\end{array}\right.$$

(44)

where ${\dot{X}}_{\mathrm{com}}\left(t\right)$ is the velocity of SIT after compensation, ${\dot{X}}_{\mathrm{com}}^{+}\left(t\right)$ is the absolute value of ${\dot{X}}_{\mathrm{com}}\left(t\right)$ when $\dot{X}\left(t\right)$ takes the maximum value ${\dot{X}}^{+}\left(t\right)$ within each vibration cycle.

When the delay compensation method for hybrid simulator achieves good force compensation effect, $\left|e_{\mathrm{F}}\right|$ and $\left|e_{\mathrm{V}}\right|$ should tend towards 0, and $e_{\mathrm{R}}$ should be close to 1.

4.1. Modelling of Switching Compensator with Variable Gain

As time of the ground hybrid simulation system increases, the differences between velocity and recovery force of hybrid simulator and the ideal system become larger without delay compensation, and the system energy will continue to increase. If fixed gain is used for compensation (FGC for short), there will be better tracking performance in the early simulation stage and gradually increasing differences in the later. Although the feedback adjustment ability can vary with time conducting linearized gain for compensation (named as LGC), the effectiveness of this compensation method is also limited when there is a significant nonlinear difference between force and velocity. In terms of the above analysis, a switching compensator with variable gain (SCVG) is constructed, which combines a linearized gain model and a variable gain model based on the energy observer, the latter of which can dynamically adjust with the change of system energy.

${\epsilon }_{\mathrm{h}}\left(k\right)$ reflects the dynamic energy fluctuations during the movement of SITS, which can be drawn as

$${\epsilon }_{\mathrm{h}}\left(k\right)=E\left(k\right)X\left(k\right)-{F^{\prime }}_x\left(k\right)X_{\mathrm{R}}\left(k\right)$$

(45)

Because the continuous oscillation of ${\epsilon }_{\mathrm{h}}\left(k\right)$ during the movement of SITS, its vibration frequency is twice the natural frequency of SIT. If it is directly introduced into differential negative feedback, it is not conducive to system stability. Therefore, integrating ${\epsilon }_{\mathrm{h}}\left(k\right)$ to obtain the feedback gain, the relationship between $F_{\mathrm{h}}\left(k\right)$ and ${\epsilon }_{\mathrm{h}}\left(k\right)$ is

$$\begin{array}{l}{K}_{\mathrm{h}}\left(k\right)=K_{\epsilon 0}+f_{\epsilon }\\ f_{\epsilon }=\left\{\begin{array}{l}{K}_{\epsilon 1}t , t\le t_{\epsilon }\\ K_{\epsilon 2}{\displaystyle \sum _k\left|{\epsilon }_{\mathrm{h}}\left(k\right)\right|}, t>t_{\epsilon }\end{array}\right.\end{array}$$

(46)

where $K_{\epsilon 0}$ is static gain, $K_{\epsilon 1}$ is variable gain, $t_{\epsilon }$ is switching time.

According to Equation (37), it can be seen that when $t\le t_{\epsilon }$, $K_{\mathrm{h}}\left(k\right)$ varies linearly and $K_{\mathrm{h}}\left(k\right)$ also changes continuously with the change of system energy when, which is related to and not entirely dependent on ${\epsilon }_{\mathrm{h}}\left(k\right)$.

When the hybrid simulation system converges, there are

$${\epsilon }_{\mathrm{h}}\left(\infty \right)=\underset{k\to \infty }{\lim }{\epsilon }_{\mathrm{h}}\left(k\right)=\underset{t\to \infty }{\lim }{\epsilon }_{\mathrm{h}}\left(t\right)=0$$

(47)

In this circumstances, $t_{\epsilon }$ is finite and ${\displaystyle \sum _k\left|{\epsilon }_{\mathrm{h}}\left(k\right)\right|}$ is bounded, so there exists $\left(C_{\mathsf{\epsilon }\mathrm{l}},C_{\mathsf{\epsilon }\mathrm{h}}\right)\in ℝ$ that satisfies

$$\begin{array}{l}\min \left\{K_{\mathsf{\epsilon }0},K_{\mathsf{\epsilon }0}+K_{\mathsf{\epsilon }1}{t}_{\mathsf{\epsilon }},K_{\mathsf{\epsilon }0}+K_{\mathsf{\epsilon }1}{\displaystyle \sum _k\left|{\epsilon }_{\mathrm{h}}\left(k\right)\right|}\right\}=C_{\mathsf{\epsilon }\mathrm{l}}\le K_{\mathrm{h}}\left(k\right)\le C_{\mathsf{\epsilon }\mathrm{h}}\\ =\max \left\{K_{\mathsf{\epsilon }0},K_{\mathsf{\epsilon }0}+K_{\mathsf{\epsilon }1}{t}_{\mathsf{\epsilon }},K_{\mathsf{\epsilon }0}+K_{\mathsf{\epsilon }1}{\displaystyle \sum _k\left|{\epsilon }_{\mathrm{h}}\left(k\right)\right|}\right\}\end{array}$$

(48)

For the negative feedback imported by SCVG, $C_{\mathsf{\epsilon }\mathrm{l}}$ is replaced by $K_{\mathrm{h}}$ in Equation (39). In addition, $\max \left({\displaystyle \sum _k\left|{\epsilon }_{\mathrm{h}}\left(k\right)\right|}\right)$ can be acquired through iterative observation, which in turn provides constraints for $K_{\epsilon 0}$ and $K_{\epsilon 1}$, which are also stability constraints for parameters tuning in the discrete domain.

4.2. Velocity Estimation Based on Tracking Differentiator

In the delay compensation approach shown in Fi.6, it is necessary to retrieve the displacement differential signal of hybrid simulator. If the position signal is directly differentiated in discrete time domain, significant noise will be involved and cause distortion of the velocity information, thereby leading to poor compensation performance. Owing to the capacity of optimal tracking differentiator in discrete-time form to rapidly track input signals without overshoot, it has advantages in eliminating signal vibration and improving differentiation quality [33,34,35,36]. Therefore, the discrete optimal tracking differentiator is applied to estimate the SITS velocity.

The actual displacement signal $X_{\mathrm{R}}\left(k\right)$ are decomposed into $x_{\mathrm{r}1}\left(k\right)$ and $x_{\mathrm{r}2}\left(k\right)$ after passing through the tracking differentiator, where $x_{\mathrm{r}1}\left(k\right)$ tracks $X_{\mathrm{R}}\left(k\right)$ and $x_{\mathrm{r}2}\left(k\right)={\dot{x}}_{\mathrm{r}1}\left(k\right)$, thus $x_{\mathrm{r}2}\left(k\right)={\widehat{\dot{X}}}_{\mathrm{R}}\left(k\right)$. The discrete optimal tracking differentiator is given by

$$\left\{\begin{array}{l}u=f_{td}\left(X_{\mathrm{R}}\left(k\right),x_{\mathrm{r}1}\left(k\right),x_{\mathrm{r}2}\left(k\right),r,h\right)\\ x_{\mathrm{r}1}\left(k+1\right)=x_{\mathrm{r}1}\left(k\right)+x_{\mathrm{r}2}\left(k\right)\cdot T\\ x_{\mathrm{r}2}\left(k+1\right)=x_{\mathrm{r}2}\left(k\right)+u\cdot T\end{array}\right.$$

(49)

where $r$ is velocity factor, which determines the tracking speed of $X_{\mathrm{R}}\left(k\right)$. The larger r, the shorter the system step response time, and the faster the command tracking speed. $h$ is filtering factor, which suppresses the noise in $X_{\mathrm{R}}\left(k\right)$. $f_{td}\left(X_{\mathrm{R}}\left(k\right),x_{\mathrm{r}1}\left(k\right),x_{\mathrm{r}2}\left(k\right),r,h\right)$ is the Optimal control synthesis function, expressed as follows

$$\left\{\begin{array}{l}d=rh\\ d_0=dh\\ y=x_{\mathrm{r}1}\left(k\right)-X_{\mathrm{R}}\left(k\right)+h\cdot x_{\mathrm{r}2}\left(k\right)\\ a_0=\sqrt{d^2+8r\left|y\right|}\\ a=\left\{\begin{array}{l}{x}_{\mathrm{r}2}\left(k\right)+\frac{a_0-d}{2}sign\left(y\right),\left|y\right|>d_0\\ x_{\mathrm{r}2}\left(k\right)+\frac{y}{h}, \left|y\right|\le d_0\end{array}\right.\\ u=\left\{\begin{array}{l}-r\cdot sign\left(a\right), \left|a\right|>d\\ -\frac{ra}{d}, \left|a\right|\le d\end{array}\right.\end{array}\right.$$

(50)

5. Simulation Study

In order to verify the correctness of the theoretical analysis above and the effectiveness of the proposed delay compensation method, simulation studies are carried out sequentially. Firstly, the comparative analyses are conducted on the SITS stability of 1D undamped vibration under different delays. Then, the response characteristic of SITS identified by the ARX model is demonstrated. Finally, the effectiveness of different delay compensation means on hybrid simulator of ground capture system for SIT is compared. Parameters required for simulation study are shown in Table 1.

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

Parameter	Value	Unit
$V_0$	0.05	m/s
$m_{\mathrm{e}}$	500	Kg
$K_{\mathrm{tc}}$	20000	N/m
${\tau }_{\mathrm{s}}$	0.002	s
$T_{\mathrm{c}}$	0.004	s
$B_{\mathrm{tc}}$	50	N·s/m
${\overline{A}}_{\mathrm{R}}$	0.025	m
${\overline{\omega }}_{\mathrm{R}}$	6π	rad/s
$T$	0.002	s
$N_{\mathrm{R}}$	1501

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

Parameter	Value	Unit
$V_0$	0.05	m/s
$m_{\mathrm{e}}$	500	Kg
$K_{\mathrm{tc}}$	20000	N/m
${\tau }_{\mathrm{s}}$	0.002	s
$T_{\mathrm{c}}$	0.004	s
$B_{\mathrm{tc}}$	50	N·s/m
${\overline{A}}_{\mathrm{R}}$	0.025	m
${\overline{\omega }}_{\mathrm{R}}$	6π	rad/s
$T$	0.002	s
$N_{\mathrm{R}}$	1501

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

(1) A hybrid simulator is proposed based on SCVG after analyzing the instability of simulator in current SPS systems caused by three types of delays. The stability constraints of the proposed hybrid simulator are obtained adopting the Routh criterion, and the simulation study is employed to verify the effectiveness and stability of the proposed compensation method.

(2) The proposed SCVG consists of variable gain and velocity estimator, which modifies input force of dynamic model through switching control. Velocity estimation is achieved by approximating tracking of differential signal of the actual displacement applying the optimal tracking differentiator. The nonlinear gain model automatically adjusts the feedback gain according to the current observed energy variation of the system.

(3) The peak errors of dynamic input, velocity, and recovery coefficient adopting SCVG in the last cycle before the end of simulation are 0.29N, 0.028mm/s, and 0.0005, respectively. Compared with FGC and LGC, the performance improvements are at least 50.00%, 82.50%, and 97.54%.

(4) The proposed SCVG demonstrates excellent compensation effects for the 1D motion of the hybrid simulator in simulation of ground capture test systems. In future work, we aim to extend the proposed compensator to 6-DoF spatial motion, further enhancing the stability of the hybrid simulator. The future research is to develop a universal hybrid simulator with high-fidelity in the ground capture mission trials.