A Novel Closed-Loop Structure for Drag-Free Control Systems with ESKF and LQR

Abstract

Space-borne gravitational wave detection satellite confronts many uncertain perturbations, such as solar pressure, dilute atmospheric drag, etc. To realize an ultra-static and ultra-stable inertial benchmark achieved by a test-mass (TM) being free to move inside a spacecraft (S/C), the drag-free control system of S/C requires super high steady-state accuracies and dynamic performances. The Active Disturbance Rejection Control (ADRC) technique has a certain capability in solving problems with common perturbations, while there is still room for optimization in dealing with the complicated drag-free control problem. When faced with complex noises, the steady-state accuracy of the traditional control method is not good enough and the convergence speed of regulating process is not fast enough. In this paper, the optimized Active Disturbance Rejection Control technique is applied. With the extended state Kalman filter (ESKF) estimating the states and disturbances in real time, a novel closed-loop control structure is designed by combining the linear quadratic regulator (LQR) and ESKF, which can satisfy the design targets competently. The comparative analysis and simulation results show that the LQR controller designed in this paper has a faster response and a higher accuracy compared with the traditional nonlinear state error feedback (NSEF), which uses a deformation of weighting components of classical PID. The new drag-free control structure proposed in the paper can be used in future gravitational wave detection satellites.

1. Introduction

Gravitational wave astronomy provides a new tool to explore black holes, dark matter, the early universe, and the evolution of the universe. To detect gravitational waves in space, a strategy involves deploying multiple satellites in mega-satellite formations to measure tiny changes in the relative distances between satellites when the gravitational wave goes through. However, the challenge lies in the weak characteristics (to the order of ${10}^{-21}$) in the changes caused by gravitational waves, which can be easily influenced by extraneous perturbations and noise. To address this issue, the typical approach is to employ drag-free satellites: an inner test mass is shielded and free-falls along the geodesic of spacetime, and the outer satellite counteracts non-conservative forces and tracks the test mass in a sensitive axis. This setup creates an ultra-static and ultra-stable platform, with the test mass serving as an inertial reference for the measurement of relative distances in space. Using this approach, it is possible to obtain accurate and reliable measurements of distance variation. And the variation represents the effect of gravitational waves.

In the drag-free control loop, there are many factors that can impact the effect of the controller, such as external environmental disturbances, sensor measurement noise, process noise, and other inevitable disturbances and noises. Models of these disturbances and noises are hard to be built precisely, which makes it difficult to determine the appropriate models for compensation. To address this challenge, the Extended State Observer (ESO) method has emerged as an effective solution for modeling, estimating, and identifying disturbances in the drag-free control loop. Using the ESO method, it is possible to obtain estimates of the disturbances, which can help the design of the controller. This method has significant implications for the development of advanced control systems for drag-free satellites, thereby improving performance and reliability.

The Active Disturbance Rejection Control (ADRC) technique, proposed by Han [1,2], combines the “anti-disturbance” and “model independence” of PID control with the idea of the state observer. The Extended State Observer (ESO) is the core of ADRC, providing a way to estimate and compensate for disturbances and uncertainties. Huang [3,4] demonstrated the design method and proof of convergence for nonlinear ESO of second- and third-order systems, showing that it can achieve fast convergence without oscillation, even in the presence of model uncertainty and disturbances. However, the complexity of nonlinear ESO increases with the growth of the number of parameters, making tuning more challenging. Despite this drawback, the effectiveness of nonlinear ESO in mitigating disturbances and uncertainties makes it a promising technique for advanced control systems in various applications. Gao [5] proposed a parameter design method for linearized ESOs based on bandwidth, which effectively reduces the design threshold and improves the convenience of the application. Yang [6] analyzed the observation error of the ESO for different forms of disturbances and concludes that when the disturbance is bounded or its derivative is bounded, ESO can effectively estimate it, and the observation error is bounded. In addition, Jin [7], Chen [8], and Gan [9] analyzed the stability of the ESO using different methods, and Shao [10] analyzed high-order ESOs by adding higher-order derivative of disturbance as the extended state. Although increasing the extended order can effectively reduce the estimation error of each state, increasing the order and bandwidth simultaneously also affects the high-frequency noise suppression effect. Therefore, a trade-off between the expansion order and bandwidth is needed to balance the estimation accuracy and the high-frequency noise suppression effect. The Extended State Kalman filter (ESKF) proposed by Xue [11] combines the advantages of both extended state observer and Kalman filter to filter the noise and estimate the system state and disturbance in dealing with nonlinear systems with strong nonlinearity, large initial estimation error, and severe noise. The Extended State Kalman filter provides a potential solution to the problem of disturbance identification for the drag-free control of gravitational wave detection satellites, when the conventional filtering methods are not sufficient to estimate the disturbance.

The Linear Quadratic Regulator (LQR) is a widely used engineering tool in the aerospace industry due to its ability to achieve optimal control under specific performance requirements. Its simplicity in design has made it a popular choice for various applications, such as quadrotors [12,13,14,15,16], hypersonic vehicles [17,18], airborne remote gimbal [19], and satellite formation problems [20,21] etc.

In the field of control engineering, the Linear Quadratic Regulator (LQR) is widely used for linear problems. To ensure the effectiveness of the controller, an accurate linear model must be established or a nonlinear model should be linearized prior to the LQR design. In cases where low control performance requirements are sufficient, typically disturbances and noise are not handled directly but are instead compensated for through control. In order to improve the robustness and disturbance resistance of the traditional LQR controller, Lu [12] introduced the Extended State Observer (ESO) to estimate random low-frequency disturbances and the estimation of ESO is used by LQR. Attitude control of spacecraft with low precision requirements, disturbances, and noise are often not preprocessed and are only compensated for through LQR controllers. For systems with higher performance requirements, such as the six-degree-of-freedom attitude control system described in Ref. [17], a combination of ESO and LQR is used to achieve higher control accuracy and stronger disturbance resistance compared to using LQR control alone. In Lin’s research [19], a standard nonlinear ESO was employed in combination with LQR to estimate and compensate for multi-source perturbations, ultimately improving the control of LQR for uncertain systems. While these studies successfully combined an extended state observer with LQR control and achieved some improvement, they focused solely on the estimation and compensation of perturbations without considering the suppression of noise.

The accuracy required for drag-free and attitude control in gravitational wave detection is crucial, so the impact of noise must be considered. Previous methods are insufficient in dealing with the noise affecting gravitational wave detection satellites, and cannot estimate the perturbations effectively. Those methods also suffer from longer setting times. To design a successful control system, it is necessary to develop effective strategies to reduce the noise impact on control performance. In this study, we propose a novel approach that combines ESKF with LQR control. We use the state and disturbance estimated by ESKF as the input information for the controller, ensuring optimized control. Our analysis and simulations show that this new approach outperforms traditional solutions. It effectively shortens the adjustment time, reduces the number of oscillations, compensates for disturbances, and suppresses noise, ultimately achieving the desired design specifications.

The paper is organized as follows: in Section 2, a dynamic model of a single test mass drag-free satellite is established. In Section 3 and Section 4, the design process and calculation methods of the ESKF and LQR are presented, respectively. In Section 5, the performance of the system using the LQR controller and NSEF controller is compared through numerical simulations, indicating that the overall performance of the LQR controller is superior to that of the NSEF controller when using ESKF as the estimation method. The conclusion is given in Section 6.

2. Dynamics Modeling

This paper takes a single test mass, a drag-free satellite in geocentric orbit as the research subject, as is shown in Figure 1, where C indicates the center of mass of an object, h represents the sensitive cavity, which is fixed to the satellite, then the position vector ${\mathbf{r}}_h$ from the center of the sensitive cavity to the center of mass (CoM) of the satellite is constant, and the position vector of test mass relative to the satellite is $\mathbf{r}={\mathbf{r}}_h+{\mathbf{r}}_{rel}$. The relative translation equations of motion in the inertial system are first transformed into the satellite body coordinate system, similarly the relative attitude equations of motion are projected into the TM body coordinate system, which is illustrated in Figure 2. Then a comprehensive drag-free satellite dynamics model can be established as follows.

$${\ddot{\phi }}_{\mathrm{sc}}={\mathrm{I}}_{\mathrm{sc}}^{-1}\left({\mathbf{T}}_{\mathrm{Csc}}+w_{T_{\mathrm{Csc}}}+{\mathbf{T}}_{Dsc}\right),$$

(1)

$$\begin{array}{cc}\hfill {\ddot{\mathbf{r}}}_{rel}^h& =\frac{1}{m_{tm}}\left({\mathbf{F}}_{Gtm}^h+{\mathbf{F}}_{Dtm}^h+{\mathbf{F}}_{SCtm}^h\right)-\frac{1}{m_{sc}}\left({\mathbf{F}}_{Gsc}^h+{\mathbf{F}}_{Csc}^h+{\mathbf{F}}_{Dsc}^h+{\mathbf{F}}_{TMsc}^h\right)\hfill \\ \hfill & -2{\mathit{\omega }}_{sc}^h\times {\dot{\mathbf{r}}}_{rel}^h-{\mathit{\omega }}_{sc}^h\times \left({\mathit{\omega }}_{sc}^h\times \left({\mathbf{r}}_h^h+{\mathbf{r}}_{rel}^h\right)\right)-{\dot{\mathit{\omega }}}_{sc}^h\times \left({\mathbf{r}}_h^h+{\mathbf{r}}_{rel}^h\right),\hfill \end{array}$$

(2)

$$\begin{array}{cc}\hfill {\ddot{\phi }}_{rel}^{tm}& ={\left({\mathbf{I}}_{tm}\right)}^{-1}\left[-\left({\mathit{\omega }}_{rel}^{tm}{\mathit{\omega }}_{sc}^{tm}\right)\times \left({\mathbf{I}}_{tm}\left({\mathit{\omega }}_{rel}^{tm}+{\mathit{\omega }}_{sc}^{tm}\right)\right)\right]+{\left({\mathbf{I}}_{tm}\right)}^{-1}\left[{\mathbf{T}}_{Gtm}^{tm}+{\mathbf{T}}_{Dtm}^{tm}+{\mathbf{T}}_{SCtm}^{tm}\right]\hfill \\ \hfill & -A_{TS}{\dot{\mathit{\omega }}}_{sc}^{sc}-A_{TS}{\mathit{\omega }}_{sc}^{sc}\times {\mathit{\omega }}_{rel}^{tm}\hfill \end{array}$$

(3)

where $w_{T_{\mathrm{Csc}}}$ indicates input noise, here assume that all input forces and moments acting on the satellite, as well as the measurement output of the sensor, are subject to noise, ${\mathbf{T}}_{Dsc}$ denotes the disturbance moment to the satellite, ${\mathbf{F}}_{Dtm}$ and ${\mathbf{F}}_{Dsc}$ denote the test mass and the disturbance force on the satellite, respectively, and ${\mathbf{T}}_{Dtm}$ indicates the disturbance moment to the test mass. The $sc$, $tm$, $rel$, C, D, and G subscripts indicate the S/C, Test Mass, measurements RELated to the sensitive cavity, Control command, Disturbance, and Gravity. The h superscript indicates the components in the sensitive cavity frame, and $A_{TS}$ is the coordinate transformation from the satellite frame to the test mass frame.

The dynamic model is expressed in the form of state space equations which are presented as follows:

$$\begin{array}{cc}\hfill {\dot{X}}_0& =A_0{X}_0+B_0\left(u+w+f\right)\hfill \\ \hfill Y& =C_0{X}_0+d\hfill \end{array}$$

(4)

where

$$\begin{array}{c}{A}_0=\left[\begin{array}{cccccc}{0}_3& I_3& 0_3& 0_3& 0_3& 0_3\\ 0_3& 0_3& 0_3& 0_3& 0_3& 0_3\\ 0_3& 0_3& 0_3& I_3& 0_3& 0_3\\ 0_3& 0_3& -\frac{K_{trans}}{m_{tm}}& -\frac{D_{trans}}{m_{tm}}& 0_3& 0_3\\ 0_3& 0_3& 0_3& 0_3& 0_3& 0_3\\ 0_3& 0_3& 0_3& 0_3& I_{tm}^{-1}{K}_{rot}& I_{tm}^{-1}{D}_{rot}\end{array}\right],B_0=\left[\begin{array}{ccc}{0}_3& 0_3& 0_3\\ I_{sc}^{-1}& 0_3& 0_3\\ 0_3& 0_3& 0_3\\ 0_3& -\frac{I_3}{m_{sc}}& 0_3\\ 0_3& 0_3& 0_3\\ -I_{sc}^{-1}& 0_3& I_{tm}^{-1}\end{array}\right],\hfill \\ C_0=\left[\begin{array}{cccccc}{I}_3& 0_3& 0_3& 0_3& 0_3& 0_3\\ 0_3& 0_3& I_3& 0_3& 0_3& 0_3\\ 0_3& 0_3& 0_3& 0_3& I_3& 0_3\end{array}\right]\hfill \end{array}$$

where $K_{trans}$, $D_{trans}$, $K_{rot}$, $D_{rot}$ are the coupling coefficient matrices for translation and rotation, respectively. u is the system control variable, w is input noise, d is measurement noise, and f represents the total perturbation affecting the system, including the known part and the unmodeled part.

The Gravitational Wave Detector-TianQin requires detection satellites in deep-space orbit. The main disturbance on the satellite comes from solar pressure. To ensure a steady power supply and minimize fluctuations in the satellite’s internal thermal environment for ultra-stability, the drag-free satellite uses a body-attached battery array.

The expression of this perturbation is shown below:

$${\overrightarrow{F}}_r=-k{C}_R{\rho }_{SR}\left(\frac{S_R}{m}\right){\overrightarrow{r}}_s$$

(5)

where $C_R$ indicates surface reflection coefficient, normally 1–1.44, ${\rho }_{SR}$ indicates the solar pressure near the Earth, $4.56\times {10}^{-6}\mathrm{N}/{\mathrm{m}}^2$, $\left(\frac{S_R}{m}\right)$ is the surface-to-mass ratio of the spacecraft, $S_R$ is the projected area of the spacecraft facing the sun, ${\overrightarrow{r}}_s$ is the unit vector indicating the direction from the center of the Earth to the sun. The sun exposure factor, denoted as k, is assumed to be 1 for the light area and 0 for the ground shadow area. The amplitude spectral density of the solar pressure on the satellite is shown in Figure 3.

In a drag-free satellite with a single test mass, the displacement measurements between the CoM of the TM and the CoM of the satellite, as well as the attitude measurements of the TM relative to the satellite, are obtained by an inertial sensor. The attitude of the satellite is determined through a star sensor, and the micro-propulsion provides the necessary control forces and moments to maintain the desired position and attitude of the satellite. At present, the typical micro-propulsion systems have a noise power spectral density of $1\times {10}^{-6}\mathrm{N}/\sqrt{\mathrm{Hz}}$ under open-loop condition [22]. Their corresponding power spectral densities are presented in Figure 4, Figure 5 and Figure 6.

Based on the data of Ref. [23], it can be inferred that electrostatic actuation noise can be treated as white noise in the frequency band needed for the detection of gravitational waves by the TianQin detector.

In the case of capacitive displacement sensors, the noise levels for displacement measurements are equal in the x, y, and z directions, while the noise levels for angle measurements in the $\theta$ direction are one order of magnitude lower than those in the $\eta$ and $\varphi$ directions [23]. Figure 6 displays the measurement noise in each direction.

The spectral density curves of perturbation and noise are given in Figure 3, Figure 4, Figure 5 and Figure 6, which will be used as the basis for the modeling and simulation calculations later.

3. Extended State Kalman Filter Design

Achieving high accuracy of relative displacement and relative attitude control in a noisy and disturbed environment requires multiple steps, including disturbance estimation, noise suppression, and state control.

The Extended State Kalman Filter (ESKF) can estimate nonlinear uncertainty. In cases of initial error, uncertain dynamics (perturbation), and bounded noise, the perturbation is estimated and compensated for by the extended state, and noise effects can be suppressed. This paper uses ESKF to estimate disturbance forces and moments, such as solar pressure on the satellite and anomalous electromagnetic forces and moments on the test mass. First, we present the design scheme of ESKF, then apply it to uncertain disturbance estimation in drag-free control.

3.1. Extended State Kalman Filter

For the following discrete system containing uncertain perturbations

$$\left\{\begin{array}{c}{W}_{k+1}=A_k{W}_k+B_{k}f\left(W_k,k\right)+w_k\hfill \\ Y_k=C_k{W}_k+n_k\hfill \end{array}\right.,k=0,1,2,\dots ,$$

(6)

where $W_k$ is the system state, $A_k,B_k$ are system matrices, $C_k$ is the measurement matrix, $f\left(W_k,k\right)$ is the nonlinear uncertain part in the system (6), $w_k,n_k$ are the process noise and measurement noise, respectively, and $Y_k$ is the system measurement output. Treating $f\left(W_k,k\right)$ as an additional state variable $f_k$, which is then estimated and compensated for. The extended system is described as

$$\left\{\begin{array}{c}\left[\begin{array}{c}{W}_{k+1}\\ f_{k+1}\end{array}\right]={A_k}^{\prime }\left[\begin{array}{c}{W}_k\\ f_k\end{array}\right]+{B_k}^{\prime }{G}_k+\left[\begin{array}{c}{w}_k\\ 0\end{array}\right]\hfill \\ Y_k={C_k}^{\prime }\left[\begin{array}{c}{W}_k\\ f_k\end{array}\right]+n_k\hfill \end{array}\right.$$

(7)

where $A_k^{\prime }=\left[\begin{array}{cc}{A}_k& B_k\\ 0& I\end{array}\right]$, $B_k^{\prime }=\left[\begin{array}{c}0\\ I\end{array}\right]$, $C_k^{\prime }=\left[\begin{array}{cc}{C}_k& 0\end{array}\right]$, $f_k=f(W_k,k)$, $G_k\stackrel{\Delta }{=}{f}_{k+1}-f_k$, assume $w_k$, $n_k$ are unrelated zero-mean Gaussian random series and $E\left(n_k{n_k}^{\mathrm{T}}\right)\le R_k$, $E\left(w_k{w_k}^{\mathrm{T}}\right)\le S_k$, $E\left(\left[\begin{array}{c}{W}_0-{\widehat{W}}_0\\ f_0-{\widehat{f}}_0\end{array}\right]{\left[\begin{array}{c}{W}_0-{\widehat{W}}_0\\ f_0-{\widehat{f}}_0\end{array}\right]}^{\mathrm{T}}\right)\le P_0$, ${\widehat{W}}_0$ is the estimation of $W_0$, ${\widehat{f}}_0$ is the initial value of the nominal part of $f\left(W_k,k\right)$, $P_0$ is a known constant matrix, and $E\left(G_i^2\right)\le {\overline{q}}_i,i=1,2,\dots$, $q_i$ is bounded.

According to the classical state observer theory, the extended state observer for the extended state Equation (7) is shown below, where $X_{k+1}=\left[\begin{array}{c}{W}_{k+1}\\ f_{k+1}\end{array}\right]$

$${\widehat{X}}_{k+1}=A_k{\widehat{X}}_k+B_k{\widehat{G}}_k-K_k\left(Y_k-C_k{\widehat{X}}_k\right)$$

(8)

Based on the given initial estimation ${\widehat{X}}_0$ and initial value of covariance matrix $P_0$, we obtain the middle value of the estimated quantity ${\widehat{X}}_k^{-}$ and update the value of the covariance matrix ${P_k}^{-}$

$${\widehat{X}}_k^{-}=A_k^{\prime }{\widehat{X}}_{k-1}+B_k^{\prime }{u}_{k-1}+B_e{\widehat{G}}_k$$

(9)

$$\begin{array}{cc}\hfill {P_k}^{-}& =(1+\theta )A_k^{\prime }{P}_{k-1}{A_k^{\prime }}^{\mathrm{T}}+\left(1+\frac{1}{\theta }\right)Q_{1,k-1}+Q_{2,k-1}\hfill \end{array}$$

(10)

where $\theta =\sqrt{\frac{\mathrm{tr}\left(Q_{1,0}\right)}{\mathrm{tr}\left(P_0\right)}}$ is used to decouple the cross terms of estimation error and uncertainty, $Q_{1,0}=4B_e{Q}_0{B_e}^{\mathrm{T}}$, $Q_{2,0}=B_k^{\prime }{S}_k{B_k^{\prime }}^{\mathrm{T}}$, $S_k$ is the variance of $w_k$, $R_k$ is the variance of $n_k$, ${\widehat{G}}_k$ is the estimation of $G_k$, whose value is calculated by ${\widehat{G}}_k=sat({\overline{G}}_k,\sqrt{q_i})$, where

$$sat\left(f,b\right)=\{\begin{array}{cc} b& f>b\\ \begin{array}{c} f\\ -b\end{array}& \begin{array}{c}b>f>-b\hfill \\ -b>f\hfill \end{array}\end{array}$$

(11)

by calculating $K_k$, update the estimation ${\widehat{X}}_k$ and the covariance matrix $P_k$

$$K_k={P_k}^{-}{C_k^{\prime }}^{\mathrm{T}}{\left({C^{\prime }}_k{P_k}^{-}{{C^{\prime }}_k}^{\mathrm{T}}+R_k\right)}^{-1}$$

(12)

$${\widehat{X}}_k={\widehat{X}}_k^{-}+K_k\left(Y_k-{C^{\prime }}_k{\widehat{X}}_k^{-}\right)$$

(13)

$$P_k=(I-K_k{C}_k^{\prime }){P_k}^{-}{(I-K_k{C}_k^{\prime })}^{\mathrm{T}}+K_k{R}_k{K_k}^{\mathrm{T}}$$

(14)

After that, calculate the control variable $u_0$ based on the error between the state estimate and the reference, then the calculated control variable $u_0$ and the estimated value $\widehat{f}$ of the disturbance are used to calculate the final control variable u, and the flow chart is shown in Figure 7.

3.2. Extended State Design of Drag-Free Control System

The complex space environment presents a challenge in accurately modeling and describing perturbations affecting on satellite and test masses. Conventional control methods are model-dependent, and their performance can be severely affected by the inaccuracies of the perturbation model. To address this issue, Section 3.1 proposes a method to expand the perturbations into new states and create a new filter model. Then the perturbation can be compensated and the noise can be suppressed.

The uncertain disturbance term is $f={[T_{Dsc},F_{Dsc},T_{Dtm}]}^{\mathrm{T}}$. $T_{Dsc}$ represents the disturbance moments that the satellite is subjected to, $F_{Dsc}$ represents the disturbance forces that the satellite is subjected to, and $T_{Dtm}$ represents the disturbance moment that TM is subjected to. Those are treated as extended states. The state vector is taken as $X={[{\phi }_{sc},{\dot{\phi }}_{sc},r_{rel},{\dot{r}}_{rel},{\phi }_{rel},{\dot{\phi }}_{rel},T_{Dsc},F_{Dsc},T_{Dtm}]}^{\mathrm{T}}$, $Y={[{\phi }_{sc},r_{rel},{\phi }_{rel}]}^{\mathrm{T}}$ indicates related attitude and displacement that can be measured.

According to Section 3.1, the extended state differential equation is obtained as

$$\begin{array}{c}\dot{X}=AX+B\left(u+w\right)+B_e\dot{f}\hfill \\ Y=CX+d\hfill \end{array}$$

(15)

where $A=\left[\begin{array}{cc}{A}_0& B_0\\ 0& I\end{array}\right],B=\left[\begin{array}{c}{B}_0\\ 0\end{array}\right],C=\left[\begin{array}{cc}{C}_0& 0\end{array}\right]$, $B_e={[0,0,1]}^{\mathrm{T}}$, symbols has the same physical meaning as in Equation (4).

The extended state Equation (15) is discretized to obtain the discrete difference equation model (16)

$$\begin{array}{cc}\hfill X_{k+1}& =A_d{X}_k+B_d\left(u+w\right)+B_d{G}_k\hfill \\ \hfill Y_k& =C_d{X}_k+d_k\hfill \end{array}$$

(16)

where $G_k\stackrel{\Delta }{=}{f}_{k+1}-f_k$.

The ESKF is designed according to the flowchart shown in Figure 7 to estimate the kinematic state parameters and the unknown disturbances, while implementing feedback control. Where $q_i={\left(\max \left|f_{i+1}-f_i\right|\right)}^2$, $Q_0=i\mathrm{diag}\left(q_i\right)$, i is the number of state vector.

4. Feedback Controller Design

Once the estimation of states and perturbations has been successfully obtained, the next step is the selection of an appropriate controller. The Linear Quadratic Regulator (LQR) is a control strategy designed to minimize a cost function. The optimal control law is obtained through the design of the state feedback controller, which allows for the completion of closed-loop optimal control in a fast, stable, and accurate manner.

The performance index for LQR control reflects the requirements for both state and control quantities, and the cost function used in this paper is

$$J=\sum _0^n\left(x{\left(k\right)}^{\mathrm{T}}Qx\left(k\right)+u{\left(k\right)}^{\mathrm{T}}Ru\left(k\right)\right)$$

(17)

The weighting matrix Q is semi-positive definite and R is positive definite, which are set as a diagonal matrix in the subsequent simulation. For the first term in the cost function J, each component is required to be small in the control process. The larger the weight in Q means the stricter the constraint on the components; while the second term in the cost function indicates the requirement for the control output, which is weighted according to the different characteristics of each component.

The Ricatti equation $PA+A^{\mathrm{T}}P-PB{R}^{-1}{B}^{\mathrm{T}}P+Q=0$ is used to calculate P, then based on $K=R^{-1}{B}^{\mathrm{T}}P$, the feedback gain matrix is calculated, the control law of LQR is chosen as $u\left(k\right)=-K\widehat{x}\left(k\right)$.

In order to maximize the effective time of gravitational wave detection and minimize the output of thrusters (actuators), it is necessary to ensure that the detection satellite is maintained in a free and stable flight for as long a time as possible. To achieve this goal, the adjustment transition process of drag-free control needs to be as fast as possible, allowing the detection system to quickly reach an ultra-static and ultra-stable state. Our research shows that the LQR control strategy can satisfy the optimal control requirements, enabling a rapid transition to a steady state and achieving the required control accuracy. By appropriately designing the weights of each component, in addition to the values of Q and R, it is possible to achieve a more refined control.

Nonlinear error feedback is used in classical Active Disturbance Rejection Control [1], mainly by rewriting the weighting of classical PID control into a nonlinear combination, as shown in Equation (18)

$$fal\left(e,\alpha ,\delta \right)=\{\begin{array}{c} \frac{e}{{\delta }^{\alpha -1}},\left|e\right|\le \delta \hfill \\ \mathrm{sgn}\left(e\right){\left|e\right|}^{\alpha },\left|e\right|>\delta \hfill \end{array}$$

(18)

It is a continuous power series function with linear segments near the origin, $\delta >0$ denotes the length of the interval of the linear segment, and e indicates the amount of error.

This function $fal\left(\right)$ is characterized by increasing the gain when the error is small and using a small gain when the error is large, which prevents high-frequency chattering due to excessive gain calculated when the error is small [24,25]. The control law for nonlinear error feedback is

$$\begin{array}{cc}\hfill u& ={\gamma }_{1}fal\left(e_1,{\alpha }_1,\delta \right)+{\gamma }_{2}fal\left(e_2,{\alpha }_2,\delta \right)+{\gamma }_{3}fal\left(e_3,{\alpha }_3,\delta \right)\hfill \end{array}$$

(19)

where $e_1\left(k\right)=r-\widehat{x}\left(k\right),e_2\left(k\right)={\displaystyle \sum _0^k}\left(r-\widehat{x}\left(k\right)\right),e_3=\dot{r}-\widehat{\dot{x}}\left(k\right)$. Parameter triples ${\gamma }_1,{\gamma }_2,{\gamma }_3$ determine the final control variable, the value of $\delta$ is generally selected as the sampling time for discrete systems, the value of $\alpha$ satisfies $\alpha \in \left(0,1\right)$.

A simulation program using the combination of ESKF and LQR is developed, then the comparison results of the combination of ESKF+NSEF and the combination of ESKF+LQR is presented and analyzed in the next section.

5. Simulation Analysis

According to the discrete model (16), a block diagram of the drag-free satellite control system is designed and presented in Figure 8. The measurement mechanism provides information on the attitude angle of the satellite, the displacement, and the attitude of the test mass relative to the center of the cavity. The output command of actuators and output measurement information of the sensors are both inputs of ESKF. The ESKF estimates the attitude, angular velocity, displacement, velocity, disturbance forces, and disturbance moments. By selecting an appropriate controller with the ESKF, high-accuracy anti-disturbance control of the drag-free satellite/test mass dynamic system is achieved.

First, we consider the LQR controller. The ESKF filter design discussed in Section 3 was simulated in MATLAB with a time step of 0.01 s. The ESKF consists of 27 states, which are divided into three groups of nine states each. After the ESKF output was stable and tracking accurately. The LQR control algorithm was then employed for designing the control law, and the corresponding simulation results were obtained.

Initial conditions of simulation:

The perturbation forces and moments of the satellite are modeled as constant superposition sinusoidal perturbations with phase differences in each axial direction. Specifically: ${\mathbf{F}}_{Dsc}=\left[\begin{array}{c}-12.8+7.7\times sin\left({\omega }_{d}t\right)\hfill \\ -12.8+7.7\times sin\left({\omega }_{d}t+\frac{2\pi }{3}\right)\hfill \\ -12.8+7.7\times sin\left({\omega }_{d}t+\frac{4\pi }{3}\right)\hfill \end{array}\right]\times {10}^{-7}\left(\mathrm{N}\right)$, the disturbance moment to the satellite is modeled as ${\mathbf{T}}_{Dsc}=\left[\begin{array}{c}-1.2+6.6\times sin\left({\omega }_{d}t\right)\hfill \\ -1.2+6.6\times sin\left({\omega }_{d}t+\frac{2\pi }{3}\right)\hfill \\ -1.2+6.6\times sin\left({\omega }_{d}t+\frac{4\pi }{3}\right)\hfill \end{array}\right]\times {10}^{-6}\left(\mathrm{N}·\mathrm{m}\right)$, moment of disturbance to the test mass is ${\mathbf{T}}_{D\mathrm{tm}}=\left[\begin{array}{c}-1.2+7.7\times sin\left({\omega }_{d}t\right)\hfill \\ -1.2+7.7\times sin\left({\omega }_{d}t+\frac{2\pi }{3}\right)\hfill \\ -1.2+7.7\times sin\left({\omega }_{d}t+\frac{4\pi }{3}\right)\hfill \end{array}\right]\times {10}^{-12}\left(\mathrm{N}·\mathrm{m}\right)$, where ${\omega }_d=1.2\times {10}^{-3}\mathrm{Hz}$. The simulation program will achieve real-time estimation of the extended state for the above perturbations, and the results of the error analysis of the estimation are given later. The expectation of the input noise of the thrusters providing force and moment to the satellite are $1\times {10}^{-9}\mathrm{N}/\sqrt{\mathrm{Hz}}$ and $1\times {10}^{-9}\mathrm{N}·\mathrm{m}/\sqrt{\mathrm{Hz}}$, respectively. The input noise expectation of the electrostatic actuator providing the test mass control torque is $1\times {10}^{-15}\mathrm{N}·\mathrm{m}/\sqrt{\mathrm{Hz}}$, the expectation of the measurement noise of the star sensor providing satellite attitude measurement is set as $1\times {10}^{-7}\mathrm{rad}/\sqrt{\mathrm{Hz}}$. And the expectation of displacement measurement noise of inertial sensor is $1\times {10}^{-9}\mathrm{m}/\sqrt{\mathrm{Hz}}$, the expectation of attitude measurement noise is $1\times {10}^{-8}\mathrm{rad}/\sqrt{\mathrm{Hz}}$. To meet the requirements of gravitational wave detection, the control loop’s design objective is set as follows: the setting time should be less than 1 min, and the amplitude spectral density of the relative displacement between test mass and satellite should both be less than ${10}^{-8}\mathrm{m}/\sqrt{\mathrm{Hz}}$ within the detection frequency band.

The performance of the combination of LQR and ESKF is evaluated initially. Subsequently, the performance of using the combination of LQR and ESKF is compared with the performance of the combination of NSEF and ESKF in terms of control accuracy and setting time. The intrinsic reasons for any differences observed are analyzed, and recommendations for engineering design are provided.

5.1. ESFK+LQR

By utilizing the ESKF-estimated states as input, we designed an LQR controller based on Equation (17). The initial value of state estimation was ${\widehat{X}}_0={[0,0,0,0,0,0,0,0,0]}^{\mathrm{T}}$, and the values of Q and R were selected as $R=1\times {10}^{-4}\mathrm{diag}(I_3,10I_3,I_3)$. It should be noted that a larger value of Q can facilitate faster convergence of the states for the LQR controller.

5.1.1. Satellite Attitude

Figure 9 demonstrates the effectiveness of the ESKF and LQR controllers in controlling the attitude angle of the satellite. The results show that the attitude angle was successfully controlled from the initial ${\phi }_{sc}={[7\times {10}^{-5},0,0]}^{\mathrm{T}}\mathrm{rad}$ to $\pm 3\times {10}^{-8}\mathrm{rad}$ for all three attitude angles with a setting time of about 7 s, achieving the desired control target. The estimation error of the attitude angle was measured to be $\pm 2\times {10}^{-8}\mathrm{rad}$, while the estimation error of the attitude angular velocity and disturbance moment were $\pm 2\times {10}^{-10}\mathrm{rad}/\mathrm{s}$ and $\pm 1\times {10}^{-5}\mathrm{N}·\mathrm{m}$, respectively.

The amplitude spectral density curves for each attitude angle control error are presented in Figure 10. Based on the results shown in Figure 10, it can be observed that within the measurement bandwidth, the amplitude spectral density of the satellite’s attitude control error conforms to the design requirements of ${10}^{-7}\mathrm{rad}/\sqrt{\mathrm{Hz}}$.

5.1.2. Test Mass and Satellite Relative Displacement

The results in Figure 11 show that the relative displacement was successfully controlled from $r_{rel}={[0.0005,0.0009,-0.0006]}^{\mathrm{T}}\mathrm{m}$ to $\pm 2\times {10}^{-9}\mathrm{m}$ for all three axes within an adjustment time of 20 s, achieving the desired control target. The position estimation error was measured to be $\pm 2.2\times {10}^{-9}\mathrm{m}$, while the velocity estimation error and disturbance estimation error were $\pm 4\times {10}^{-10}\mathrm{m}/\mathrm{s}$ and $\pm 4\times {10}^{-8}\mathrm{N}$, respectively.

The amplitude spectral density curves for each axial direction are presented in Figure 12. Based on the observations from Figure 12, it can be concluded that the kinematic indexes of the translation within the frequency band satisfy the design requirements and achieve ${10}^{-8}\mathrm{m}/\sqrt{\mathrm{Hz}}$.

5.1.3. Relative Attitude between Test Mass and Satellite

Figure 13 shows the simulation results of ESKF and LQR controller dealing with the relative attitude between test mass and satellite. Based on the results from Figure 13, it can be concluded that the ESKF and LQR controllers designed in this paper were capable of controlling the relative attitude between the test mass and the satellite from ${\phi }_{rel}={[0,4\times {10}^{-5},0]}^{\mathrm{T}}\mathrm{rad}$ to $\pm 5\times {10}^{-8}\mathrm{rad}$ for all three axes. The control result was found to be essentially oscillation-free, with an adjustment time of about 6s, achieving the desired control target. The estimation error of the attitude angle was measured to be $\pm 5\times {10}^{-8}\mathrm{rad}$, while the estimation errors of the attitude angular velocity and disturbance moment were $\pm 6\times {10}^{-8}\mathrm{rad}/\mathrm{s}$ and $\pm 2\times {10}^{-11}\mathrm{N}·\mathrm{m}$, respectively.

The amplitude spectral density curves for each attitude angle are presented in Figure 14. The amplitude spectral densities are lower than $5\times {10}^{-8}\mathrm{rad}/\sqrt{\mathrm{Hz}}$ in the whole frequency band.

5.2. Comparison Evaluation

To demonstrate the efficacy of the designed method, we present a simulation to compare ESKF+LQR with ESKF+NSEF. This allows us to illustrate the functions of our objectives and evaluate their performance against an established method. Simulation results is given in Figure 15, Figure 16 and Figure 17

Since the convergence facility of the $fal\left(\right)$ function used in NSEF is mainly based on the value of $\alpha$, while the proof of certain physical relation is extremely complicated, the value of $\alpha$ is determined by rule of thumb. Based on the method in Ref. [26], we take ${\alpha }_1={\alpha }_2={\alpha }_3=0.5$ in simulation. The simulation results of relative displacement of the x-axis are shown in Figure 16, while the control accuracy is satisfied, and the convergence process of ESKF+LQR was found to be much faster than that of ESKF+NSEF.

The RMS of the two approaches is listed in

Table 1. RMS of two different control approaches.

Axis	ESKF+LQR	ESKF+NSEF
satellite attitude (pitch)	6.3076 × ${10}^{-9}$ rad	7.4918 × ${10}^{-9}$ rad
satellite attitude (yaw)	7.2104 × ${10}^{-9}$ rad	5.9713 × ${10}^{-9}$ rad
satellite attitude (roll)	7.4293 × ${10}^{-9}$ rad	4.7966 × ${10}^{-9}$ rad
relative displacement (x)	4.5918 × ${10}^{-10}$ m	8.9468 × ${10}^{-10}$ m
relative displacement (y)	4.5422 × ${10}^{-10}$ m	9.3517 × ${10}^{-10}$ m
relative displacement (z)	4.7442 × ${10}^{-10}$ m	9.4479 × ${10}^{-10}$ m
relative attitude (pitch)	1.0859 × ${10}^{-10}$ rad	2.7201 × ${10}^{-8}$ rad
relative attitude (yaw)	1.0480 × ${10}^{-10}$ rad	3.1758 × ${10}^{-8}$ rad
relative attitude (roll)	1.0299 × ${10}^{-10}$ rad	3.2608 × ${10}^{-8}$ rad

. Additionally, the LQR controller was observed to produce negligible oscillations in the system, which can be advantageous in terms of reducing energy consumption and extending the effective time of gravitational wave signal detection.

6. Conclusions

Integrating perturbations into states vector, the ESKF method demonstrates the capability to accurately estimate the state and disturbance of the drag-free satellite dynamics, and effectively suppressing handling noise. This approach lays the foundation for the controller to achieve accurate adjustment of the system state, and it is recommended to consider the ESKF method as a viable alternative for estimating uncertain disturbances in future drag-free engineering designs.

The LQR controller’s feedback parameter matrix can be rigorously derived using the generalized index to ensure optimal control performance. By combining the reasonable values of the state covariance matrix, the designed gain matrix guarantees the convergence speed of the system state, while the accurate estimation of ESKF enhances the relative attitude control accuracy by an order of magnitude.

For the relative displacement, the adjustment time of NSEF fails to reach the control task requirements. This is mainly due to the relatively slow control process of the relative kinematic parameters between the satellite and the test mass in NSEF, and the lack of a complete theoretical method to adjust and optimize the nonlinear error feedback parameters.

In summary, the ESKF+LQR control approach enables the system to reach a steady state rapidly and smoothly, thereby increasing the free flight time available for gravitational wave observation. This is more in line with the desired observation duration for gravitational wave detection. In comparison to the combination of ESKF+NSEF, the combination of ESKF+LQR is capable of adjusting to the reference target swiftly and reducing oscillations, leading to a reduction in thruster mass consumption.