Attitude Stabilization of a Satellite with Large Flexible Elements Using On-Board Actuators Only

Abstract

Attitude control of a satellite with three flexible elements is considered. Control torque is developed by a set of reaction wheels, which are installed on the central hub of the satellite. The flexible elements are large, so the control torque constraints must be taken into account. In the paper, a control algorithm based on a linear-quadratic regulator is studied. The asymptotic stability of this control is shown. The choice of the control parameters is based on the closed form solution of the corresponding algebraic Riccati equation, which is supplemented by the linear matrix inequality. To increase the convergence rate, particle swarm optimization is used to tune the control parameters.

1. Introduction

Consideration of elastic deformations is crucial to ensure high pointing precision for large space structures (for example, the ISS [1]), and, especially, for satellites with antennas [2], solar panels [3] and robotic manipulators [4]. Generally, the abovementioned elements have relatively low natural frequencies and damping ratios, so elastic deformations of such elements affect the whole system dynamics and could even lead to instabilities. Therefore, high-precision attitude control of flexible spacecraft has to suppress vibrations.

Numerous control design methods have been adapted and applied for attitude control of flexible spacecrafts. First of all, this concerns classical methods, such as proportional–integral–derivative (PID) controller [3], a linear quadratic regulator (LQR) [5] and sliding mode (SM) control [6]. In [7] attitude trajectory, the tracking problem is solved by means of LMI-based gain-scheduled H-infinity control. The SM control concept related to attitude maneuvers of a flexible spacecraft was further developed. In particular, Cao et al. [8] presented a robust attitude control law based on a novel nonsingular terminal sliding surface and operating in the face of actuator uncertainty and uncertain spacecraft dynamics. An SM adaptive fault-tolerant attitude tracking control algorithm is proposed for flexible spacecraft with partial loss of actuator effectiveness, unknown inertia parameters and external disturbances [9]. Additionally, using the singular perturbation theory [10,11], the equations of attitude motion of a flexible satellite are decomposed into fast and slow subsystems, and in both cases, the SM control laws are constructed relying on the latter one. Moreover, a novel component synthesis vibration suppression method is presented [12]. The authors note that the proposed active vibration suppression technique is consistent with such controllers as PD controllers, SM controllers and model predictive controls. The latter was used to accomplish large-angle single axis rotational maneuvers of a flexible satellite [13] and attitude/spin maneuvers of a spacecraft equipped with a large rotating flexible reflector (NASA Soil Moisture Passive Active mission) [2]. A number of research papers are devoted to the important issue of ensuring the high quality of transient processes during attitude maneuvers [13,14]. In [14], a composite adaptive neural prescribed performance control ensuring the prescribed performance of transient response and attitude trajectory convergence in a preselected finite settling time is proposed. Fixed-time attitude control and stabilization using a neural network is considered in [15]. Fuzzy-logic optimal control was investigated in [16]. Thus, there are a huge variety of control approaches for flexible spacecraft attitude stabilization.

Control algorithms usually need parameter tuning. In the paper, particle swarm optimization (PSO) is used. In [17], the state-dependent Riccati equation (SDRE) attitude control technique is used within the scope of three flexible satellites formation flying, and the PSO algorithm tuned the pulse-width and pulse-frequency (PWPF) modulator, generating on–off thruster commands. The paper in [18] presented two combined LQR-PSO and SM-PSO attitude control methods, within which the PSO algorithm optimized the variable angle between the rigid hub and the payload. Note that not only the control law parameters appear as the subject of optimization. For instance, the kinematics of a rigid satellite undergoing constrained slew maneuvers using reaction wheels are tuned by an inverse-dynamics PSO approach [19].

Besides derivation of the attitude control law, there is the problem of obtaining its input information, such as attitude, angular velocity and the values of flexible deformations. The available set of actuators and sensors as well as their location obviously affects the choice of attitude control and vibration suppression strategies for a satellite with flexible appendages. In the last two decades, to solve the latter of these aforementioned problems, scientists have paid great attention to the use of smart structures, especially piezoelectric materials, as actuators or sensors [3,10,11,20,21,22]. For instance, Song and Agrawal [23] utilize PWPF modulation for producing the on–off thruster firing sequence required for attitude maneuver execution. At the same time, smart sensors and smart actuators with positive position feedback (PPF) control are used to actively suppress vibrations at the flexible appendage. However, the problem of simultaneous angular motion control of the satellite and vibration suppression in flexible elements, by means of actuators and sensors located only on the central body (hub), is still relevant. In this case, the clusters of momentum-exchanging devices are often employed as actuators, such as a system of four pyramid-mounted control moment gyros [1], ortho-skew construction of the four reaction wheels [24], and systems of six reaction wheels [12]. Concerning measuring devices and sensors for the motion estimation of spacecraft with flexible elements, Ivanov et al. [25] provided a relevant description as well as an overview of methods for the vibration determination. In [26], the comparison between Extended Kalman filter and Unscented Kalman filter is made through a Monte-Carlo simulation. In both cases, magnetometer and sun sensor measurements are used.

The attitude control in this paper is based on two assumptions that come from practical applications: the actuators and sensors are located on the central rigid body, and the data on the eigenmodes and eigen frequencies are rather inaccurate. Moreover, usually the performance of the on-board computer is relatively low, so the state vector cannot include a large number of variables. The principal idea is rather trivial—the control utilizes the rigid model of the spacecraft and ignores the flexibility. Thus, the control algorithm becomes as simple as possible for the on-board, almost real-time, implementation. It utilizes the estimation of the satellite position, attitude and angular velocity only. However, this simplicity creates a set of problems that be solved. First is the stability problem. The system model includes a large number of degrees of freedom, and the control effects directly only a few of them; moreover, most are completely ignored during the control calculation stage. So, a stability analysis of the system is needed, and in the paper, the theoretical basis of this approach is presented, and sufficient conditions are derived. Secondly, every algorithm has control parameters that allow one to tune the system behavior. Additionally, in the paper, the methodology for selection of the control parameters using PSO is provided.

In summary, this paper is devoted to the derivation of the attitude control law for a satellite with flexible appendages enabling low-frequency vibration suppression. After the problem statement (Section 2), in Section 3 a mathematical model of a large flexible satellite, consisting of a satellite hub and three flexible elements—two solar arrays and an antenna—is presented. Flexible deformations of each element are described by the corresponding eigenmodes of vibration [27,28,29]. Section 4 contains linearized equations of satellite attitude motion. Section 5 and Section 6 are dedicated to control synthesis methodology. It is assumed that actuators and sensors are located only on the hub, so the proposed control methodology does not require additional actuators on flexible elements. Also, it does not contain any information about modal amplitudes to avoid complicating the state estimation process and increasing the computational complexity of the control algorithm with their identification. Hence, a reduced (rigid) model of the spacecraft is used to obtain an attitude control law. In this context the following methods are tested: an inertia-free nonlinear attitude control algorithm derived by Sanyal et al. [30] and implemented by Posani et al. [24]; LQR [5] as well as SDRE techniques; the SDRE algorithm tuned by an input-shaping technique to reduce undesired elastic oscillations [31]. In the present paper, the LQR-PSO strategy is presented. The LQR is a well-known approach. However, in the present paper, the control is based on the reduced model (for a rigid body with fewer degrees of freedom) while being applied to the full one, which includes deformation modes. Direct implementation of the LQR in this case may lead to instability, so proper selection of the control parameters must be achieved. In Section 5, the sufficient conditions for asymptotic stability are derived, and an explicit solution of the algebraic Riccati equation (ARE) for this system is found. In Section 6, PSO searches for the optimal values of the LQR parameters minimizing the system’s degree of stability. To improve the calculation effectiveness, the results of Section 5 are used. In Section 7, an illustration is presented. The appendices contain the d’Alembert principle for the whole system application and general forces calculation.

2. Problem Statement

The problem of attitude control of a satellite with three flexible elements (FE) is considered. The spacecraft is a rigid central hub with two flexible panels and an antenna cantilever attached to it (Figure 1).

The actuators (reaction wheels) and sensors (star tracker and angular velocity sensor) are installed on the central hub only. Therefore, there is neither direct action on nor direct measurements of the flexible deformations. The goal of the control system is to stabilize the hub in the inertial space and damp oscillations in the FEs.

3. Equations of Motion

There are two models, which are used in the paper: the full nonlinear model for the numerical modelling and linear model for the control algorithm synthesis, which is obtained from the first one. This section presents the nonlinear model.

The following reference frames are used:

• $OXYZ$ is the nonrotating frame; its origin coincides with Earth’s center of mass, $OZ$ is perpendicular to the equatorial plane, $OX$ is directed to the vernal equinox point corresponding to a given epoch (e.g., J2000);
• $O_{s}xyz$ is the body-fixed frame; its origin lies in the satellite hub center of mass ($S$), and its axes coincide with its principal axes of inertia;
• $O_k{x}_k{y}_k{z}_k,\hspace{0.17em}k=\overline{1,3}$ are the flexible-element fixed frames with origin in the center of mass of the corresponding undeformed flexible element; axes are the principal axes of inertia of the undeformed flexible element.

The points of the satellite hub $m_{s,i}$, solar panel $m_{p,i}$ and antenna $m_{a,i}$ positions (in Figure 1) are defined as follows:

$$\begin{array}{l}{R}_{s,i}=R_s+r_{s,i},\hspace{0.17em}\\ R_{k,i}=R_k+r_{k,i}+{\rho }_{k,i},\hspace{0.17em}k=\overline{1,3},\end{array}$$

(1)

where $R_s,\hspace{0.17em}\hspace{0.17em}{R}_k$ are the radius vectors of $S$ and $F{E}_k$ centers of mass, respectively, given in $OXYZ$, $r_{s,i},\hspace{0.17em}\hspace{0.17em}{r}_{k,i}$ are the radius vectors of $S$ and $F{E}_k$ i-th points with respect to $O_{s}xyz$, $O_k{x}_k{y}_k{z}_k$, respectively, and ${\rho }_{k,i}$ is the displacement of $F{E}_k$ i-th point due to deformation, respectively. It is assumed that ${\rho }_{k,i}\ll r_{k,i}\ll R_k,\hspace{0.17em}\hspace{0.17em}k=\overline{1,3}$.

Normal modes are used for deformation definition. Point displacements due to deformations for each $F{E}_k$ are [27,28,29]

$${\rho }_{k,i}=A_{k,i}{q}_k,$$

(2)

where $q_k\left(t\right)={\left(q_1\left(t\right),\hspace{0.17em}\hspace{0.17em}\dots \hspace{0.17em}\hspace{0.17em},q_{n_k}\left(t\right)\right)}^T$ are the vectors of modal coordinates, $n_k$ is the number of modes taken into account, $A_{k,i}$ is the $3\times n_k$ matrix of the mode shapes. The j-th column $\left(1\le j\le n_k\right)$ of these matrices defines the flexible displacements of the i-th point of the $F{E}_k$ (panel or antenna) caused by its j-th normal mode.

There are various approaches to derive the equations of motion of a flexible multibody system [27,32,33]. In this paper, the nonlinear model was developed using the d’Alembert [34] principle for each part of the satellite. The corresponding equations for the hull are [29,35]

$$S_s\left(\begin{array}{c}{\ddot{R}}_s\\ {\dot{\omega }}_s\end{array}\right)=N_s$$

(3)

where

$$S_s=\left(\begin{array}{cc}{m}_s{I}_{3\times 3}& 0\\ 0& J_s\end{array}\right),\hspace{0.17em}\hspace{0.17em}{N}_s=\left(\begin{array}{c}{F}_s\\ M_s-{\omega }_s\times J_s{\omega }_s\end{array}\right).$$

(4)

Here, ${\omega }_s$ is the absolute angular velocity of the hub, $F_s$ and $M_s$ are the net force and torques acting on the hub (including reaction forces in the joint points), respectively, and $m_s$ and $J_s$ are the mass and the inertia tensor of the hub, respectively. $I_{n\times n}$ is the n-by-n identity matrix.

The corresponding equations for the $F{E}_k$ are [29,35,36]

$$S_k\left(\begin{array}{c}{\ddot{R}}_k\\ {\dot{\omega }}_k\\ {\ddot{q}}_k\end{array}\right)=N_k,$$

(5)

where

$$S_k=\left(\begin{array}{ccc}{m}_k{I}_{3\times 3}& -{\left[m_k{A}_k{q}_k\right]}_{\times }& m_k{A}_k\\ {\left[m_k{A}_k{q}_k\right]}_{\times }& {\tilde{J}}_k& {\displaystyle \sum _{i\in F{E}_k}{\left[m_{k,i}\left(r_{k,i}+{\rho }_{k,i}\right)\right]}_{\times }{A}_{k,i}}\\ m_k{A}_k^T& -{\displaystyle \sum _{i\in F{E}_k}{A}_{k,i}^T{\left[m_{k,i}\left(r_{k,i}+{\rho }_{k,i}\right)\right]}_{\times }}& I_{n_k\times n_k}\end{array}\right),$$

(6)

$$N_k=\left(\begin{array}{c}{F}_k-m_k{\omega }_k\times {\omega }_k\times A_k{q}_k-2m_k{\omega }_k\times A_k{\dot{q}}_k\\ {\displaystyle \sum _{i\in F{E}_k}\left(r_{k,i}+{\rho }_{k,i}\right)\times F_{k,i}}-{\omega }_k\times {\tilde{J}}_k{\omega }_k-2{\displaystyle \sum _{i\in F{E}_k}{m}_{k,i}\left(r_{k,i}+{\rho }_{k,i}\right)\times {\omega }_k\times {\dot{\rho }}_{k,i}}\\ -{\Omega }_{n_k\times n_k}{q}_k+{\displaystyle \sum _{i\in F{E}_k}{m}_{k,i}{A}_{k,i}^T{F}_{k,i}}-{\displaystyle \sum _{i\in F{E}_k}{m}_{k,i}{A}_{k,i}^T{\omega }_k\times {\omega }_k\times \left(r_{k,i}+{\rho }_{k,i}\right)}-2{\displaystyle \sum _{i\in F{E}_k}{m}_{k,i}{A}_{k,i}^T{\omega }_k\times {\dot{\rho }}_{k,i}}\end{array}\right).$$

(7)

Here, ${\omega }_k$ is the absolute angular velocity of the $F{E}_k$, $F_k$ and $M_k$ are the net force and torques acting on the $F{E}_k$ (including reaction forces in the joint points), respectively, and $m_k$ and ${\tilde{J}}_k$ are, respectively, the mass and the inertia tensor of the $F{E}_k$. The last one accounts for the deformations, also.

$$A_k={\displaystyle \sum _{i\in F{E}_k}{m}_{k,i}{A}_{k,i}}.$$

(8)

$${\Omega }_{n_k\times n_k}=diag\left(\begin{array}{cccc}{\mathsf{\Omega }}_{k,1}^2& {\mathsf{\Omega }}_{k,2}^2& \dots & {\mathsf{\Omega }}_{k,n_k}^2\end{array}\right)$$

(9)

${\mathsf{\Omega }}_{k,i}$ is the i-th eigen frequency of the $F{E}_k$. The notation ${\left[r\right]}_{\times }$ for a vector $r={\left(\begin{array}{ccc}{r}_1& r_2& r_3\end{array}\right)}^T$ is

$${\left[r\right]}_{\times }=\left(\begin{array}{ccc}0& -r_3& r_2\\ r_3& 0& -r_1\\ -r_2& r_1& 0\end{array}\right).$$

Finally, the equation for the full state vector accelerations $\ddot{x}={\left(\begin{array}{ccccc}{\ddot{R}}_s^T& {\dot{\omega }}_s^T& {\ddot{q}}_1^T& {\ddot{q}}_2^T& {\ddot{q}}_3^T\end{array}\right)}^T$ (the details can be found in the Appendix A) is

$$\tilde{S}\ddot{x}=\tilde{N},$$

(10)

where the total dynamics matrix $\tilde{S}$ and right part vector $\tilde{N}$ are

$$\tilde{S}=\left(\begin{array}{cccc}{S}_s+{\displaystyle \sum _{k=1}^3{W}_{k,1}^T{S}_k{W}_{k,1}}& W_{11}^T{S}_1{W}_{12}& W_{21}^T{S}_2{W}_{22}& W_{31}^T{S}_3{W}_{32}\\ W_{12}^T{S}_1{W}_{21}& W_{12}^T{S}_1{W}_{12}& 0& 0\\ W_{22}^T{S}_2{W}_{21}& 0& W_{22}^T{S}_2{W}_{22}& 0\\ W_{32}^T{S}_3{W}_{31}& 0& 0& W_{32}^T{S}_3{W}_{32}\end{array}\right),$$

(11)

$$\tilde{N}=\left(\begin{array}{c}{N}_s+{\displaystyle \sum _{k=1}^3\left(W_{k,1}^T{N}_k-W_{k,1}^T{S}_k{T}_k\right)}\\ W_{12}^T{N}_1-W_{12}^T{S}_1{T}_1\\ W_{22}^T{N}_2-W_{22}^T{S}_2{T}_2\\ W_{32}^T{N}_3-W_{32}^T{S}_3{T}_3\end{array}\right).$$

(12)

Here,

$$W_{k,1}=\left(\begin{array}{cc}{E}_{3\times 3}& -{\left[s_k+f_k\right]}_{\times }\\ 0_{3\times 3}& I_{3\times 3}\\ 0_{n_k\times 3}& 0_{n_k\times 3}\end{array}\right),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{W}_{k,2}=\left(\begin{array}{c}{0}_{3\times n_k}\\ 0_{3\times n_k}\\ E_{n_k\times n_k}\end{array}\right),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{T}_k=\left(\begin{array}{c}{\omega }_s\times {\omega }_s\times \left(s_k+f_k\right)\\ 0_{1\times 3}\\ 0_{1\times n_k}\end{array}\right).$$

(13)

Equation (10) is supplemented by the kinematic relations, which have the form

$$\begin{array}{l}{\dot{R}}_s={\mathbf{V}}_s,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\\ {\dot{\lambda }}_{0s}=-\frac{1}{2}\left({\omega }_s,\hspace{0.17em}{\lambda }_s\right),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\\ {\dot{\lambda }}_s=\frac{1}{2}\left({\lambda }_{0s}{\omega }_s+{\lambda }_s\times {\omega }_s\right),\\ {\dot{\mathbf{q}}}_k={\mathbf{V}}_{q_k},\hspace{0.17em}\hspace{0.17em}k=\overline{1,3},\end{array}$$

(14)

where ${\left(\begin{array}{cc}{\lambda }_{0s}& {\lambda }_s^T\end{array}\right)}^T$ is the attitude quaternion of the hub, ${\mathbf{V}}_{q_k}$ is the change rate vector of normal mode amplitudes for corresponding flexible elements $F{E}_k$. Equations (10) and (14) completely determine the motion of the considered system.

In comparison with [29], the current equations of motion are written relative to the center of mass of the rigid hub. In this case, if the number of flexible elements attached to the hub needs to be changed, the system (10) can be easily modified.

4. Linearized Mathematical Model

To derive the attitude control and study the asymptotic stability of the desired equilibrium, the linearized equations of satellite motion relative to its center of mass are used. Here, the required position

$$\begin{array}{l}{\mathsf{\omega }}_s={\mathbf{0}}_{3\times 1},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\mathsf{\lambda }}_s={\mathbf{0}}_{3\times 1},\\ {\mathbf{q}}_k={\mathbf{0}}_{n_k\times 1},{\mathbf{V}}_{{\mathbf{q}}_k}={\mathbf{0}}_{n_k\times 1},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}k=\overline{1,3},\end{array}$$

(15)

is considered. It is also supposed that there are no external torques and forces except control torque $u$ produced by the set of actuators installed on the hub. The reasons for this assumption will be discussed in Section 5. Also, the orbital motion must be excluded. So, the system takes the following form:

$$\left(\begin{array}{cccc}{\mathbf{S}}_{\omega }& {\mathbf{S}}_{\omega q_1}& {\mathbf{S}}_{\omega q_2}& {\mathbf{S}}_{\omega q_3}\\ {\mathbf{S}}_{\omega q_1}^T& {\mathbf{S}}_{q_1}& {\mathbf{S}}_{q_1{q}_2}& {\mathbf{S}}_{q_1{q}_3}\\ {\mathbf{S}}_{\omega q_2}^T& {\mathbf{S}}_{q_1{q}_2}^T& {\mathbf{S}}_{q_2}& {\mathbf{S}}_{q_2{q}_3}\\ {\mathbf{S}}_{\omega q_3}^T& {\mathbf{S}}_{q_1{q}_3}^T& {\mathbf{S}}_{q_2{q}_3}^T& {\mathbf{S}}_{q_3}\end{array}\right)\left(\begin{array}{c}{\dot{\mathsf{\omega }}}_s\\ {\ddot{\mathsf{q}}}_1\\ {\ddot{\mathsf{q}}}_2\\ {\ddot{\mathsf{q}}}_3\end{array}\right)=\left(\begin{array}{c}u\\ -{\Omega }_1{q}_1\\ -{\Omega }_2{q}_2\\ -{\Omega }_3{q}_3\end{array}\right),$$

(16)

$$\begin{array}{l}{\mathbf{S}}_{\omega }=\mathbf{J}={\mathbf{J}}_s+{\displaystyle \sum _{k=1}^3\left({\mathbf{J}}_k-m_k{\left[s_k+f_k\right]}_{\times }{\left[s_k+f_k\right]}_{\times }\right)}+\frac{1}{m}{\displaystyle \sum _{k=1}^3\left(m_k{\left[s_k+f_k\right]}_{\times }\right)}{\displaystyle \sum _{k=1}^3\left(m_k{\left[s_k+f_k\right]}_{\times }\right)},\\ {\mathbf{S}}_{\omega q_k}={\displaystyle \sum _{i\in F{E}_k}{m}_{k,i}{\left[r_{k,i}\right]}_{\times }{A}_{k,i}}+m_k\left({\left[s_k+f_k\right]}_{\times }-\frac{1}{m}{\displaystyle \sum _{k=1}^3{m}_k{\left[s_k+f_k\right]}_{\times }}\right)A_k,\hspace{0.17em}\hspace{0.17em}k=\overline{1,3},\\ {\mathbf{S}}_{q_k}=E_{n_k\times n_k}-\frac{m_k^2}{m}{\mathbf{A}}_k^T{\mathbf{A}}_k,\hspace{0.17em}\hspace{0.17em}k=\overline{1,3},\\ {\mathbf{S}}_{q_k{q}_l}=-\frac{m_k{m}_l}{m}{\mathbf{A}}_k^T{\mathbf{A}}_l,\hspace{0.17em}\hspace{0.17em}k\ne l=\overline{1,3}.\end{array}$$

(17)

Grouping flexible variables as $q={\left(\begin{array}{ccc}{q}_1^T& q_2^T& q_3^T\end{array}\right)}^T$, $V_q={\left(\begin{array}{ccc}{V}_{q_1}^T& V_{q_2}^T& V_{q_3}^T\end{array}\right)}^T$ and $\Omega =diag\left(\begin{array}{ccc}{\Omega }_1& {\Omega }_2& {\Omega }_3\end{array}\right)$, and taking into account the kinematic Equation (14), gives

$$\left(\begin{array}{cc}\mathbf{J}& S_{\omega q}\\ S_{\omega q}^T& S_q\end{array}\right)\left(\begin{array}{c}{\dot{\omega }}_s\\ {\dot{V}}_q\end{array}\right)=-\left(\begin{array}{cc}0& 0\\ 0& \Omega \end{array}\right)\left(\begin{array}{c}0\\ q\end{array}\right)+\left(\begin{array}{c}u\\ 0\end{array}\right),$$

(18)

where

$$S_{\omega q}=\left(\begin{array}{ccc}{S}_{\omega q_1}& S_{\omega q_2}& S_{\omega q_3}\end{array}\right),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{S}_q=\left(\begin{array}{ccc}{S}_{q_1}& S_{q_1{q}_2}& S_{q_1{q}_3}\\ S_{q_1{q}_2}^T& S_{q_2}& S_{q_2{q}_3}\\ S_{q_1{q}_3}^T& S_{q_2{q}_3}^T& S_{q_3}\end{array}\right),$$

(19)

is the matrix of natural vibration frequencies of flexible elements.

Having solved the system (18) with respect to higher derivatives and using kinematics (14), the linear equations of angular motion are [37]

$$\dot{x}=Ax+Bu,$$

(20)

where $x={\left(\begin{array}{cccc}{\omega }_s^T& V_q^T& {\lambda }_s^T& {\mathbf{q}}^T\end{array}\right)}^T$ is the state vector of SC with flexible elements,

$$\begin{array}{l}A=\left(\begin{array}{cccc}{0}_{3\times 3}& 0_{3\times n_{\Sigma }}& 0_{3\times 3}& {\mathbf{J}}^{-1}{S}_{\omega q}{\left(S_q-S_{\omega q}^T{\mathbf{J}}^{-1}{S}_{\omega q}\right)}^{-1}\Omega \\ 0_{n_{\Sigma }\times 3}& 0_{n_{\Sigma }\times n_{\Sigma }}& 0_{n_{\Sigma }\times 3}& -{\left(S_q-S_{\omega q}^T{\mathbf{J}}^{-1}{S}_{\omega q}\right)}^{-1}\Omega \\ \frac{1}{2}{E}_{3\times 3}& 0_{3\times n_{\Sigma }}& 0_{3\times 3}& 0_{3\times n_{\Sigma }}\\ 0_{n_{\Sigma }\times 3}& E_{n_{\Sigma }\times n_{\Sigma }}& 0_{n_{\Sigma }\times 3}& 0_{n_{\Sigma }\times n_{\Sigma }}\end{array}\right),\\ B=\left(\begin{array}{c}{\mathbf{J}}^{-1}\left(E_{3x3}+S_{\omega q}{\left(S_q-S_q^T{\mathbf{J}}^{-1}{S}_{\omega q}\right)}^{-1}{S}_{\omega q}^T{\mathbf{J}}^{-1}\right)\\ -{\left(S_q-S_{\omega q}^T{\mathbf{J}}^{-1}{S}_{\omega q}\right)}^{-1}{S}_{\omega q}^T{\mathbf{J}}^{-1}\\ 0_{3\times 3}\\ 0_{n_{\Sigma }\times 3}\end{array}\right).\end{array}$$

(21)

Here, the number $n_{\Sigma }={\displaystyle {\sum }_{n=1}^3{n}_k}$ denotes the sum of all modes considered in the satellite’s model. The latter system is used for the stabilization part of the control.

5. Control Synthesis

Since the control must stabilize in the inertial space, there will be two main external torques that affect the angular motion in the required equilibrium: gravity gradient and solar pressure. Both values are small (the first is about ${10}^{-3}\hspace{0.17em}\mathrm{N}\cdot \mathrm{m}$, the latter—${10}^{-4}\hspace{0.17em}\mathrm{N}\cdot \mathrm{m}$) with respect to the possible control torque (maximum is $1\hspace{0.17em}\hspace{0.17em}\mathrm{N}\cdot \mathrm{m}$). The control consists of two parts: stabilizing and compensating. The second part compensates the external torques only, while the first one is based on the linear-quadratic regulator (LQR) [37] and does not include these torques. This control part is based on the reduced linearized model with no external torques and not requiring determination of the amplitudes of the eigenmodes to calculate the control of interest.

5.1. Stabilizing Control

The LQR is based on a linear solid-state model of the satellite, so in (20), the state vector and the matrices take the form

$$x=x_s={\left(\begin{array}{cc}{\omega }_s^T& {\lambda }_s^T\end{array}\right)}^T,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}A=A_s=\left(\begin{array}{cc}{0}_{3\times 3}& 0_{3\times 3}\\ \frac{1}{2}{E}_{3\times 3}& 0_{3\times 3}\end{array}\right),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}B=B_s=\left(\begin{array}{c}{J}^{-1}\\ 0_{3\times 3}\end{array}\right).$$

(22)

The LQR minimizes the following cost function [37]

$$I={\displaystyle \underset{0}{\overset{\infty }{\int }}\left(x_s^{T}Q{x}_s+u^{T}Ru\right)dt}$$

(23)

with positive definite matrices $Q$, $R$ that are the parameters of the algorithm. The LQR has the form [37]

$$u=-R^{-1}{B}_s^{T}P{x}_s,$$

(24)

where $P$ is a solution to the algebraic Riccati equation [37]

$$A_s^{T}P+P{A}_s-P{B}_s{R}^{-1}{B}_s^{T}P+Q=0.$$

(25)

Since $A_s$ is the $6\times 6$ matrix, $P$ size is also $6\times 6$.

Since the pair $\left(\mathbf{A},\mathbf{B}\right)$ from (22) is controllable and $Q$, $R$ are positive definite matrices, the following are known [37]:

• Matrix $P$ from the LQR control law is the only positive definite solution of (25);
• The LQR provides the asymptotic stability for the linear system with matrices (22).

Let the matrix $Q$ have the form

$$Q=\left(\begin{array}{cc}{Q}_{11}& Q_{12}\\ Q_{12}^T& Q_{22}\end{array}\right),$$

(26)

where ${\mathbf{Q}}_{11}$ and ${\mathbf{Q}}_{22}$ are positive definite matrices (${\mathbf{Q}}_{11}>0$, ${\mathbf{Q}}_{22}>0$). A positive definite solution of (25) $P$ is convenient to represent as

$$\mathbf{P}=\left(\begin{array}{cc}{\mathbf{P}}_{11}& {\mathbf{P}}_{12}\\ {\mathbf{P}}_{12}^T& {\mathbf{P}}_{22}\end{array}\right),\hspace{0.17em}$$

(27)

where ${\mathbf{P}}_{11}={\mathbf{P}}_{11}^T$, ${\mathbf{P}}_{12}$, ${\mathbf{P}}_{22}={\mathbf{P}}_{22}^T$ are $3\times 3$ matrices. This representation is possible since $A_s$ and $B_s$ are block matrices with $3\times 3$ elements. So, Equation (25) is rewritten as

$$\frac{1}{2}\left(\left(\begin{array}{cc}{\mathbf{P}}_{12}^T& {\mathbf{P}}_{22}\\ {\mathbf{0}}_{3\times 3}& {\mathbf{0}}_{3\times 3}\end{array}\right)+\left(\begin{array}{cc}{\mathbf{P}}_{12}& {\mathbf{0}}_{3\times 3}\\ {\mathbf{P}}_{22}& {\mathbf{0}}_{3\times 3}\end{array}\right)\right)-\left(\begin{array}{cc}{\mathbf{P}}_{11}\mathbf{Z}{\mathbf{P}}_{11}& {\mathbf{P}}_{11}\mathbf{Z}{\mathbf{P}}_{12}\\ {\mathbf{P}}_{12}^T\mathbf{Z}{\mathbf{P}}_{11}& {\mathbf{P}}_{12}^T\mathbf{Z}{\mathbf{P}}_{12}\end{array}\right)+\left(\begin{array}{cc}{Q}_{11}& Q_{12}\\ Q_{12}^T& Q_{22}\end{array}\right)=0.$$

(28)

Here $\mathbf{Z}=J^{-1}{R}^{-1}{J}^{-1}$ is the positive definite matrix since the inertia tensor is $J>0$.

As a result, Equation (25) is represented as a system of three equations:

$$\frac{1}{2}\left({\mathbf{P}}_{12}^T+{\mathbf{P}}_{12}\right)-{\mathbf{P}}_{11}\mathbf{Z}{\mathbf{P}}_{11}+{\mathbf{Q}}_{11}=0,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\frac{1}{2}{\mathbf{P}}_{22}={\mathbf{P}}_{11}\mathbf{Z}{\mathbf{P}}_{12}+{\mathbf{Q}}_{12},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\mathbf{P}}_{12}^T\mathbf{Z}{\mathbf{P}}_{12}={\mathbf{Q}}_{22}.$$

(29)

The fourth equation is exactly the same as the second equation in the system (29):

$${\mathbf{P}}_{12}^T\mathbf{Z}{\mathbf{P}}_{11}+{\mathbf{Q}}_{12}^T=\frac{1}{2}{\mathbf{P}}_{22}=\frac{1}{2}{\mathbf{P}}_{22}^T={\left({\mathbf{P}}_{12}^T\mathbf{Z}{\mathbf{P}}_{11}\right)}^T={\mathbf{P}}_{11}\mathbf{Z}{\mathbf{P}}_{12}+{\mathbf{Q}}_{12}.$$

(30)

LQR control with (22) taken into account is

$$u=-R^{-1}{J}^{-1}\left(P_{11}{\omega }_s+P_{12}{\lambda }_s\right)=-K_{\omega }{\omega }_s-K_{\lambda }{\lambda }_s,$$

(31)

where $K_{\omega }=R^{-1}{J}^{-1}{P}_{11}$, $K_{\lambda }=R^{-1}{J}^{-1}{P}_{12}$. Since the control is based on the reduced system, it is necessary to select such matrices $K_{\omega }$, $K_{\lambda }$ that provide the asymptotic stability for the full system.

Consider the following Lyapunov function:

$$V=\frac{1}{2}\left(\begin{array}{cc}{\omega }_s^T& V_q^T\end{array}\right)\left(\begin{array}{cc}J& S_{\omega q}\\ S_{\omega q}^T& S_q\end{array}\right)\left(\begin{array}{c}{\omega }_s\\ V_q\end{array}\right)+{\lambda }_s^T{K}_{\lambda }{\lambda }_s+\frac{1}{2}{q}^T\Omega q.$$

(32)

To satisfy $V>0$, it is sufficient that $K_{\lambda }>0$. Its derivative due to the equations of motion is

$$\dot{V}={\omega }_s^T\left(u+K_{\lambda }{\lambda }_s\right)=-{\omega }_s^T{K}_{\omega }{\omega }_s.$$

(33)

If $K_{\omega }>0$, then $\dot{V}\le 0$, and $\dot{V}=0$ if and only if ${\omega }_s=0$. According to the Barbashin–Krasovski–LaSalle theorem [38], the equilibrium is asymptotically stable if there are no other trajectories on the set $\left\{\dot{V}=0\right\}$ besides the equilibrium. Such a set here is $\left\{{\omega }_s=0\right\}$; this implies

$$\begin{array}{l}{S}_{\omega q}{\dot{V}}_q=-K_{\lambda }{\lambda }_s,\\ S_q{\dot{V}}_q=-\Omega y.\end{array}$$

(34)

The right part of the first equation is a constant (${\lambda }_s=const$). In the second equation, the matrices $S_q$ and $\Omega$ are positive definite, so its solution can be represented as

$$y=H\left(\begin{array}{c}{C}_1\cos \left({\nu }_{1}t+{\phi }_{01}\right)\\ C_2\cos \left({\nu }_{2}t+{\phi }_{02}\right)\\ ⋮\\ C_N\cos \left({\nu }_{N}t+{\phi }_{0N}\right)\end{array}\right)=H\Theta \left(t\right),$$

(35)

where $H$ consists of vectors $h_i$ defined from the following equation:

$$\left(\Omega -{\nu }_i{S}_q\right)h_i=0.$$

(36)

From the first Equation (34), it follows that

$$S_{\omega q}H\ddot{\Theta }\left(t\right)=-K_{\lambda }{\lambda }_s.$$

(37)

This equation is not satisfied for any values of the constants $C_i\ne 0$, ${\phi }_{0i}$, if the frequencies ${\nu }_i$ are incommensurable and $rank\left(S_{\omega q}H\right)=3$. Thus, the condition of the asymptotic stability of the zero solution is obtained. In practice, both conditions are usually fulfilled, but in some configurations (e.g., two identical panels installed symmetrically), this approach can face difficulties.

Thus, in the case of satisfying the abovementioned conditions for the satellite’s dynamics, a sufficient condition for asymptotic stability is the positive definiteness of the quadratic forms $K_{\omega }$, $K_{\lambda }$. This means that $R^{-1}$, $J^{-1}$ and $P_{11}$ should commute, since these matrices are positive definite. If $P_{12}$ is positive definite, then $R^{-1}$, $J^{-1}$ and $P_{12}$ should commute, too. For three positive definite quadratic forms to commute, it is necessary and sufficient that they have a diagonal form in the same basis. Since the inertia tensor $J$ is defined, it determines the basis in which the remaining quadratic forms should have a diagonal form. In this case, $J=\mathbf{W}{J}^{diag}{\mathbf{W}}^T,$ where $\mathbf{W}$ is an orthogonal matrix and $J^{diag}=diag\left(J_1,J_2,J_3\right)$. Equation (20) with matrices (22) and, therefore, the system (29), has a unique solution $P>0$.

First, consider the case when $J$ is diagonal. If $Q$ and $R$ are diagonal and their diagonal elements $R_i>0,\hspace{0.17em}\hspace{0.17em}i=\overline{1,3}$, and $Q_i>0,\hspace{0.17em}\hspace{0.17em}i=\overline{1,6}$, then the matrix $\mathbf{Z}=J^{-1}{R}^{-1}{J}^{-1}$ is diagonal, i.e., $\mathbf{Z}=diag\left(z_1,z_2,z_3\right)$, where $z_i=1/\left(R_i{J}_i^2\right),\hspace{0.17em}\hspace{0.17em}i=\overline{1,3}$.

In this case, the solution of the system (29) consists of diagonal matrices ${\mathbf{P}}_{11}^{diag}$, ${\mathbf{P}}_{12}^{diag}$, ${\mathbf{P}}_{22}^{diag}$

$${\mathbf{P}}_{11}^{diag}=\left(\begin{array}{ccc}{p^{\prime }}_1& 0& 0\\ 0& {p^{\prime }}_2& 0\\ 0& 0& {p^{\prime }}_3\end{array}\right),\hspace{0.17em}\hspace{0.17em}{\mathbf{P}}_{12}^{diag}=\left(\begin{array}{ccc}{p}_1& 0& 0\\ 0& p_2& 0\\ 0& 0& p_3\end{array}\right),\hspace{0.17em}\hspace{0.17em}{\mathbf{P}}_{22}^{diag}=\left(\begin{array}{ccc}{p^{″}}_1& 0& 0\\ 0& {p^{″}}_2& 0\\ 0& 0& {p^{″}}_3\end{array}\right),$$

(38)

and it is positive definite, and hence the only stabilizing one. Indeed, a particular solution of the third equation from the system (29) is a matrix ${\mathbf{P}}_{12}$ with the following diagonal elements:

$$p_i=\pm \sqrt{Q_{i+3}/z_i}=\pm J_i\sqrt{R_i{Q}_{i+3}},\hspace{0.17em}\hspace{0.17em}i=\overline{1,3}.$$

(39)

Then, from the first equation, the diagonal elements of matrix ${\mathbf{P}}_{11}$ are

$$\begin{array}{l}{p^{\prime }}_i=\pm \sqrt{\frac{p_i+Q_i}{z_i}}=\pm \sqrt{\frac{\pm \sqrt{Q_{i+3}/z_i}+Q_i}{z_i}}=\\ =\pm J_i\sqrt{R_i\left(\pm J_i\sqrt{R_i{Q}_{i+3}}+Q_i\right)},\hspace{0.17em}\hspace{0.17em}i=\overline{1,3},\end{array}$$

(40)

and we obtain the diagonal elements of the matrix ${\mathbf{P}}_{22}$, substituting the coefficients $p_i$ and ${p^{\prime }}_i$ into the second equation:

$$\begin{array}{l}{p^{″}}_i=2z_i{p^{\prime }}_i{p}_i=2z_i\sqrt{\frac{Q_{i+3}}{z_i}\frac{p_i+Q_i}{z_i}}=2\sqrt{Q_{i+3}\left(p_i+Q_i\right)}=\\ =2\sqrt{q_{i+3}\left(\sqrt{Q_{i+3}/z_i}+Q_i\right)},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}i=\overline{1,3}.\end{array}$$

(41)

According to the Schur lemma [39], it is necessary that ${\mathbf{P}}_{11}>0$ and ${\mathbf{P}}_{22}>0$ in order for $\mathbf{P}>0.$ Hence, for the positive definiteness of the Riccati equation’s solution, it is required to choose solutions with a plus sign when extracting the roots. It remains to check the fulfillment of the condition ${\mathbf{P}}_{22}-{\mathbf{P}}_{12}^T{\mathbf{P}}_{11}^{-1}{\mathbf{P}}_{12}>0$ or ${p^{″}}_i-{\left(p_i\right)}^2/{p^{\prime }}_i>0,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}i=\overline{1,3}$ in the case of diagonal matrices. Taking into account (39), (40) and (41), the latter becomes $\sqrt{Q_{i+3}/z_i}+2Q_i>0,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}i=\overline{1,3}$. This inequality is always satisfied because the coefficients $z_i,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}i=\overline{1,3}$ and $Q_i,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}i=\overline{1,6}$ are positive.

Thus, based on (39), (40) and (41), a positive definite solution of the algebraic Riccati equation is

$$\mathbf{P}=\left(\begin{array}{cc}diag\left({p^{\prime }}_1,{p^{\prime }}_2,{p^{\prime }}_3\right)& diag\left(p_1,p_2,p_3\right)\\ diag\left(p_1,p_2,p_3\right)& diag\left({p^{″}}_1,{p^{″}}_2,{p^{″}}_3\right)\end{array}\right).$$

(42)

The linear quadratic control law (31) in this case has the form

$$u=-diag\left[{\left(\sqrt{\left(\sqrt{Q_{i+3}/R_i}{J}_i+Q_i/R_i\right)}\right)}_{i=1}^3\right]{\omega }_s-diag\left[{\left(\sqrt{Q_{i+3}/R_i}\right)}_{i=1}^3\right]{\lambda }_s.$$

(43)

Now, consider the case when $J$ is a non-diagonal positive definite matrix. Let the matrices $Q_{11},$ $Q_{22}$ and $R$ be diagonal in the same basis as the inertia tensor, i.e.,

$$\begin{array}{l}{Q}_1=\mathbf{W}{Q}_1^{diag}{\mathbf{W}}^T=\mathbf{W}diag\left(Q_1,Q_2,Q_3\right){\mathbf{W}}^T,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\\ {\mathbf{Q}}_2=\mathbf{W}{Q}_2^{diag}{\mathbf{W}}^T=\mathbf{W}diag\left(Q_4,Q_5,Q_6\right){\mathbf{W}}^T,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\\ R=\mathbf{W}{\mathbf{R}}^{diag}{\mathbf{W}}^T=\mathbf{W}diag\left(R_1,R_2,R_3\right){\mathbf{W}}^T,\end{array}$$

(44)

where the diagonal elements $r_i>0,\hspace{0.17em}\hspace{0.17em}i=\overline{1,3}$ and $q_i>0,\hspace{0.17em}\hspace{0.17em}i=\overline{1,6}$. The search for a stabilizing Riccati equation solution leads to the diagonal case. Particular solution of the third equation of the system (29) is sought in the form

$${\mathbf{P}}_{12}=\mathbf{W}{\mathbf{P}}_{12}^{diag}{\mathbf{W}}^T$$

(45)

where ${\mathbf{P}}_{12}^{diag}=diag\left(p_1,p_2,p_3\right)$ is a diagonal matrix with the corresponding positive coefficients (39) of the diagonal inertia tensor case. Substitution of (45) into the first equation of the system (29) leads to the solution

$${\mathbf{P}}_{11}=\mathbf{W}{\mathbf{P}}_{11}^{diag}{\mathbf{W}}^T,$$

(46)

where ${\mathbf{P}}_{11}^{diag}=diag\left({p^{\prime }}_1,{p^{\prime }}_2,{p^{\prime }}_3\right)$ is given by the expression (40). Similarly, the matrix ${\mathbf{P}}_{22}$ is calculated as

$${\mathbf{P}}_{22}=\mathbf{W}{\mathbf{P}}_{22}^{diag}{\mathbf{W}}^T,$$

(47)

where the diagonal elements of the matrix ${\mathbf{P}}_{22}^{diag}=diag\left({p^{″}}_1,{p^{″}}_2,{p^{″}}_3\right)$ are given by the expression (41). Here, we use commutability of the diagonal matrices, the orthogonality of the matrix $\mathbf{W}$ and the fact that for an orthogonal matrix $\mathbf{W}$ and an arbitrary positive definite diagonal matrix $\mathbf{\Xi }$, the matrix $\mathbf{W}\mathbf{\Xi }{\mathbf{W}}^T$ is positive definite.

Then, the solution of ARE (29) is represented as

$$\mathbf{P}=\left(\begin{array}{cc}\mathbf{W}& {\mathbf{0}}_{3\times 3}\\ {\mathbf{0}}_{3\times 3}& \mathbf{W}\end{array}\right)\left(\begin{array}{cc}{\mathbf{P}}_{11}^{diag}& {\mathbf{P}}_{12}^{diag}\\ {\mathbf{P}}_{12}^{diag}& {\mathbf{P}}_{22}^{diag}\end{array}\right)\left(\begin{array}{cc}{\mathbf{W}}^T& {\mathbf{0}}_{3\times 3}\\ {\mathbf{0}}_{3\times 3}& {\mathbf{W}}^T\end{array}\right)$$

(48)

and it is positive definite. Herewith, the linear-quadratic control law (31)

$$u=-R^{-1}{J}^{-1}\left(\mathbf{W}{\mathbf{P}}_{11}^{diag}{\mathbf{W}}^T\hspace{0.17em}{\omega }_s+\mathbf{W}{\mathbf{P}}_{12}^{diag}{\mathbf{W}}^T{\lambda }_s\right)=-K_{\omega }{\omega }_s-K_{\lambda }{\lambda }_s$$

(49)

stabilizes a linear system with matrices (22), and hence, the system with matrices (21). Moreover, since it is true for a linearized system, the same is true for the initial non-linear system also [38].

Thus, the possibility of stabilizing the system using the specified control law is shown, and the method for solving the Riccati equation for the considered problem is described.

5.2. Compensation Control

The provided stabilization control part ignores external torques. There are two torques that should be considered in the geostationary orbit for a satellite with large solar panels: solar radiation pressure and gravity gradient torque.

First, the general expression for external torques must be derived. The (12) has the form

$$N=\left(\begin{array}{c}{N}_s+{\displaystyle \sum _{k=1}^3{W}_{k,1}^T{N}_{k,}}\\ W_{12}^T{N}_1\\ W_{22}^T{N}_2\\ W_{32}^T{N}_3\end{array}\right), \hspace{0.17em}{N}_s=\left(\begin{array}{c}{F}_s\\ M_s\end{array}\right), N_k=\left(\begin{array}{c}{F}_k\\ {\displaystyle \sum _{i\in F{E}_k}\left(r_{k,i}+{\rho }_{k,i}\right)\times F_{k,i}}\\ {\displaystyle \sum _{i\in F{E}_k}{m}_{k,i}{A}_{k,i}^T{F}_{k,i}}\end{array}\right),\hspace{0.17em}\hspace{0.17em}k=\overline{1,3}.$$

(50)

Since

$$W_{k,1}^T{N}_k=\left(\begin{array}{c}{F}_k\\ {\left[s_k+f_k\right]}_{\times }{F}_k+{\displaystyle \sum _{i\in \mathit{F}{\mathit{E}}_k}\left(r_{k,i}+{\rho }_{k,i}\right)\times F_{k,i}}\end{array}\right),\hspace{0.17em}\hspace{0.17em}k=\overline{1,3},$$

(51)

$$W_{k,2}^T{N}_k={\displaystyle \sum _{i\in \mathit{F}{\mathit{E}}_k}{m}_{k,i}{A}_{k,i}^T{F}_{k,i},}\hspace{0.17em}\hspace{0.17em}k=\overline{1,3}$$

(52)

The resulting general force vector is

$$N=\left(\begin{array}{c}{M}_s+{\displaystyle \sum _{k=1}^3\left({\left[s_k+f_k\right]}_{\times }{F}_k+{\displaystyle \sum _{i\in F{E}_k}\left(r_{k,i}+{\rho }_{k,i}\right)\times F_{k,i}}\right)-\frac{1}{m}\left({\displaystyle \sum _{k=1}^3{m}_k{\left[s_k+f_k\right]}_{\times }}\right)F}\\ {\displaystyle \sum _{i\in F{E}_1}{m}_{1,i}{A}_{1,i}^T{F}_{1,i}}-\frac{1}{m}{m}_1{A}_1^{T}F\\ {\displaystyle \sum _{i\in F{E}_2}{m}_{2,i}{A}_{2,i}^T{F}_{2,i}}-\frac{1}{m}{m}_2{A}_2^{T}F\\ {\displaystyle \sum _{i\in F{E}_3}{m}_{3,i}{A}_{3,i}^T{F}_{3,i}}-\frac{1}{m}{m}_3{A}_3^{T}F\end{array}\right),$$

(53)

where

$$F=F_s+{\displaystyle \sum _{k=1}^3{F}_k}$$

(54)

is the net force acting upon the system. The first component of (53) is the net torque of all forces acting upon the whole system. The torques, which are taken into account in the compensation part of the control, are $M^{grav}$ and $M^{sun}$. The first one corresponds to the rigid body gravity gradient torque (see Appendix B) [40]

$$M^{grav}=3\frac{\mu }{R^3}R\times JR.$$

(55)

In this case, the resulting control torque is

$$M^{ctrl}=-M^{grav}-M^{sun}-u.$$

(56)

The solar radiation torque for symmetrical configuration of the panels is thought to be near zero and is considered as a perturbation due to the small difference between the panel and its mounting point. This control law is used for the spacecraft stabilization in the neighborhood of the required position.

6. Optimization Problem

LQR requires the choice of control parameters, i.e., diagonalized matrices $Q_{11}$, $Q_{22}$ and $R.$ These matrices significantly affect the maximum values of the control vector components and the quality of the algorithm. Here the degree of stability [41] is taken as a quality metric. The maximum affordable control torque value is the upper bound. For linear systems with constant coefficients whose equilibrium position is asymptotically stable, the degree of stability is the distance from the imaginary axis to the rightmost root of the characteristic equation. In fact, this is the exponent with the least damping.

A closed-loop system with control

$$u=-Kx,$$

(57)

where $K=\left(\begin{array}{cccc}{K}_{\omega }& 0_{3\times n}& K_{\lambda }& 0_{3\times n}\end{array}\right)$, is represented in the form

$$\dot{x}=\left(A-BK\right)x=A_{c}x.$$

(58)

Eigenvalues of matrix $A_c$ determine the transient process rate. In this case, it is also necessary to take into account the constraint on control

$$\max \left|Kx\right|\le u_{\max }.$$

(59)

Thus, the optimization problem is

$$\Phi =-\mathrm{Re}{\mu }_{\min }\to \max \hspace{0.17em}\hspace{0.17em}\mathrm{under} \max \left|Kx\right|\le u_{\max }.$$

(60)

Here, ${\mu }_{\min }$ is the eigenvalue with the minimum distance to the imaginary axis on the complex plane of characteristic values.

It is necessary to formalize the criterion (59), since as a rule, there is a surge effect, i.e., an increase in the required control at the beginning of transients. For this, the approach described in [42] is used. The Lyapunov function (32) is

$$V=\frac{1}{2}\left(\begin{array}{cccc}{\omega }_s^T& V_q^T& {\lambda }_s^T& q^T\end{array}\right)\left(\begin{array}{cccc}{S}_{\omega }& S_{\omega q}& 0_{3\times 3}& 0_{3\times N}\\ S_{\omega q}^T& S_q& 0_{N\times 3}& 0_{N\times N}\\ 0_{3\times 3}& 0_{3\times N}& 2K_{\lambda }& 0_{N\times N}\\ 0_{N\times 3}& 0_{N\times N}& 0_{N\times N}& \Omega \end{array}\right)\left(\begin{array}{c}{\omega }_s\\ V_q\\ {\lambda }_s\\ q\end{array}\right)=\frac{1}{2}{x}^{T}Hx.$$

(61)

At the initial moment $t=t_0$, we have $0.5\hspace{0.17em}{x}_0^{T}H{x}_0=a_0$, where $x_0=x\left(t_0\right)$. Due to the decreasing Lyapunov function, we obtain $0.5\hspace{0.17em}{x}^{T}Hx\le a_0$ or

$$\frac{1}{2a_0}{x}^{T}Hx\le 1.$$

(62)

Rewrite the condition (59) in the form

$$x^T{K}^{T}Kx\le u_{\max }^2.$$

(63)

Thus, it is necessary to guarantee (63) under the condition (62). This is true if and only if the matrix inequality

$$\frac{1}{u_{\max }^2}{K}^{T}K\le \frac{1}{2a_0}H$$

(64)

is satisfied. By Schur’s lemma, this is equivalent to the fact that the following matrix is positive definite.

$$\left(\begin{array}{cc}H& K^T\\ K& \frac{u_{\max }^2}{2a_0}{E}_{3\times 3}\end{array}\right)\ge 0.$$

(65)

So, if the linear feedback matrix $K$ satisfies (65), the constraint will not be violated, and simultaneously, due to the decreasing (61), the system will not leave the initial ball (62).

In this case, the cost function and restrictions have complex form and cannot be presented explicitly, depending on the problem parameters. To solve this problem, the evolutionary optimization method—particle swarm optimization (PSO) [43]—is implemented for the search of the set of optimal control parameters.

Let $x_p$ be a set of the control parameters. The PSO is based on the decision-making model of each swarm particle. The model describing the decision making of particles in a swarm turned out to be a simple and effective optimization method. The task of the swarm is to provide a minimum of the given cost function

$$\Phi \left(x_p\right):{ℝ}^D\to ℝ$$

(66)

on the search domain

$$\mathbb{U}=\left\{x_p\in {ℝ}^D\hspace{0.17em}\hspace{0.17em}|\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\eta }_{low,j}\le x_{p,j}\le {\eta }_{up,j},\hspace{0.17em}\hspace{0.17em}j=\overline{1,D}\right\}$$

(67)

defined by restrictions on the values of the $D$ parameters. Each particle $p=\overline{1,P}$ in each generation $i=\overline{1,G}$ has a certain position $x_p(i)$ and velocity $v_p(i)$. The position of the particle specifies a possible solution to the optimization problem.

The velocity allows deciding the direction of displacement to continue the search and consists of three components:

$$v_p\left(i\right)=c_{in}{v}_p\left(i-1\right)+c_{cog}\left(x_{p,best}\left(i\right)-x_p\left(i\right)\right)+c_{soc}\left(x_{p,local\hspace{0.17em}best}\left(i\right)-x_p\left(i\right)\right).$$

(68)

The position of each particle at the next iteration is determined based on its current position and velocity:

$$x_p\left(i+1\right)=x_p\left(i\right)+v_p\left(i\right)$$

(69)

The first term in (68) is the inertial component; it is responsible for the search continuation in the same direction. The second is the cognitive component; the particle desires to return to its own best position found earlier, $x_{p,best}$. The last one is the social component, which represents striving for the best position found in the particle vicinity, $x_{p,local\hspace{0.17em}best}$.

The value of the cost function ${\Phi }_{p,local\hspace{0.17em}best}$ corresponds to the best particle $p$ position $x_{p,best}$, ${\Phi }_{p,local\hspace{0.17em}best}$ corresponds to the best position in the vicinity of the particle $p$ $x_{p,local\hspace{0.17em}best}$:

$$\begin{array}{c}{\Phi }_{p,best}=\Phi \left(x_{p,best}\right)\\ {\Phi }_{p,local\hspace{0.17em}best}=\Phi \left(x_{p,local\hspace{0.17em}best}\right)\end{array},$$

(70)

and ${\Phi }_{best}$ is the global optimal solution found by the entire swarm in $i$ iterations.

The contribution of each velocity component is varied using appropriate weighting coefficients $c_{in},\hspace{0.17em}\hspace{0.17em}{c}_{cog},\hspace{0.17em}\hspace{0.17em}{c}_{soc}$. A large value of the inertial coefficient $c_{in}$ accelerates the exploration of search space and does not allow it to fall into local optima. The correct selection of social $c_{soc}$ and cognitive $c_{cog}$ coefficients allows each particle to first look for its best position, and then switch attention to improving the best position found among all of the particle’s neighbors. For each optimization problem, the coefficients should be selected individually, focusing on some empirical rules and selection methods, which are given, for example, in [43]. However, in any problem when selecting coefficients, it is necessary that the following relationship be satisfied

$$c_{in}>\frac{1}{2}\left(c_{soc}+c_{cog}\right)-1$$

(71)

to ensure the convergence of the PSO, which was shown in [44].

The search stop criterion is the fulfillment of the following conditions simultaneously:

• the cost function derivative is small (dimensionless parameter of cost function stagnation is ${\mathsf{\Phi }}_{stagn}$);
• all particles are falling into some neighborhood of the best position (dimensionless parameter of swarm stagnation is $S_{stagn}$).

Figure 2 shows a block diagram of the described algorithm for clarity.

As a result, PSO searches for the optimal values of the LQR parameters, taking into account the fulfillment of condition (65). Minimization of the degree of stability can be carried out with a different number of known modes.

Since matrices $Q_{11}$, $Q_{22}$ and $R$ are diagonal, there are nine parameters that should be found. However, matrices $K_{\omega }$ and $K_{\lambda }$ depend on the ratios of $Q_{11}$ and $Q_{22}$ to $R$, so one of these matrices can be fixed. Consider that $R=diag\left(\begin{array}{ccc}1& 1& 1\end{array}\right)$ in the basis of principle axes of the system inertia tensor. The other six parameters are defined by the PSO method. So $x_p$ (66) is $x_p={\left(\begin{array}{cccccc}{Q}_1& Q_2& Q_3& Q_4& Q_5& Q_6\end{array}\right)}_p$, which are presented in (44).

Since (65) depends on the initial condition, the LQR parameters will also depend on the initial state. So, two sets of PSO bounds are used to improve the convergence in the neighborhood of the equilibrium. The first domain is

$$\begin{array}{l}\left|{\omega }_i\right|\le {10}^{-3}rad/sec,\\ \left|{\lambda }_i\right|\le 0.5;\end{array}$$

(72)

the second one is

$$\begin{array}{l}\left|{\omega }_i\right|\le {10}^{-6}rad/sec,\\ \left|{\lambda }_i\right|\le 3\times {10}^{-4}.\end{array}$$

(73)

It is considered that satellites are deployed with rather small angular momentum, so the initial domain for the angular velocity is about 10 times more than the orbital angular velocity of the geostationary orbit, which is ${\omega }_{orb}=0.7\times {10}^{-5}\hspace{0.17em}rad/sec$. Such selection of the initial condition domains is based on the following practical consideration. The spacecraft with large FE usually separates from the launch vehicle with rather small initial angular velocity, and its FE are undeployed. The deployment of FE increases the inertia tensor by two orders and hence decreases the angular velocity by the same value due to the conservation of angular momentum. The final domain is about 10 times less than the ${\omega }_{orb}$. The control limit is $u_{\max }=1\hspace{0.17em}\hspace{0.17em}N\cdot m$. PSO parameters are given in

Table 1. Main parameters of the particle swarm method.

D	$D=6$
${\eta }_{low}$ and ${\eta }_{up}$	${\eta }_{low}^1=[{10}^2\cdot \left(1,1,1\right),\hspace{0.17em}\hspace{0.17em}0.1\cdot \left(1,1,1\right)]$ ${\eta }_{up}^1=[{10}^5\cdot \left(1,1,1\right),\hspace{0.17em}\hspace{0.17em}{10}^2\cdot \left(1,1,1\right)]$	${\eta }_{low}^2=[{10}^6\cdot \left(1,1,1\right),\hspace{0.17em}\hspace{0.17em}{10}^5\cdot \left(1,1,1\right)]$ ${\eta }_{up}^2=[{10}^9\cdot \left(1,1,1\right),\hspace{0.17em}\hspace{0.17em}{10}^8\cdot \left(1,1,1\right)]$
$P$	200
$G$	500
$c_{in}^{up},\hspace{0.17em}\hspace{0.17em}{c}_{in}^{low}$	0.9, 0.4
$c_{cog}^{up},\hspace{0.17em}\hspace{0.17em}{c}_{cog}^{low}$	2.05, 0
$c_{soc}^{up},\hspace{0.17em}\hspace{0.17em}{c}_{soc}^{low}$	2.05, 0
${\mathsf{\Phi }}_{stagn}$	0.001
$S_{stagn}$	0.005

.

So, the two following sets of LQR parameters are found:

$$Q_{11}^1=\left(\begin{array}{ccc}8.49\times {10}^5& -0.02& 0.85\\ -0.02& 8.18\times {10}^5& 4.32\times {10}^3\\ 0.85& 4.32\times {10}^3& 4.40\times {10}^3\end{array}\right)\hspace{0.17em}\hspace{0.17em}{(\mathrm{N}\cdot \mathrm{m}\cdot \mathrm{s})}^2,\hspace{0.17em}\hspace{0.17em}{Q}_{22}^1=\left(\begin{array}{ccc}0.45& 0& 0\\ 0& 0.43& 0.02\\ 0& 0.02& 0.1\end{array}\right)\hspace{0.17em}\hspace{0.17em}{(\mathrm{N}\cdot \mathrm{m})}^2$$

(74)

$$Q_{11}^2=\left(\begin{array}{ccc}4.34\times {10}^7& 17.9& 444\\ 17.9& 3.92\times {10}^7& 2.12\times {10}^6\\ 444& 2.12\times {10}^6& 1.18\times {10}^6\end{array}\right)\hspace{0.17em}\hspace{0.17em}{(\mathrm{N}\cdot \mathrm{m}\cdot \mathrm{s})}^2,\hspace{0.17em}\hspace{0.17em}{Q}_{22}^2=\left(\begin{array}{ccc}2.31\times {10}^5& 1.35& 0.54\\ 1.35& 1.00\times {10}^5& -4.79\times {10}^3\\ 0.54& -4.79\times {10}^3& 1.86\times {10}^5\end{array}\right)\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{(\mathrm{N}\cdot \mathrm{m})}^2$$

(75)

The corresponding feedback matrices $K_{\omega }$ and $K_{\lambda }$ are

$$K_{\omega }^1=\left(\begin{array}{ccc}412& 0& 0\\ 0& 404& 16\\ 0& 16& 120\end{array}\right)\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{N}\cdot \mathrm{m}\cdot \mathrm{s},\hspace{0.17em}\hspace{0.17em}{K}_{\lambda }^1=\left(\begin{array}{ccc}0.67& 0& 0\\ 0& 0.66& 0\\ 0& 0& 0.32\end{array}\right)\hspace{0.17em}\hspace{0.17em}\mathrm{N}\cdot \mathrm{m}$$

(76)

$$K_{\omega }^2=\left(\begin{array}{ccc}1.02& 0& 0\\ 0& 0.89& 0.03\\ 0& 0.03& 0.39\end{array}\right)\times {10}^4\hspace{0.17em}\mathrm{N}\cdot \mathrm{m}\cdot \mathrm{s},\hspace{0.17em}\hspace{0.17em}{K}_{\lambda }^2=\left(\begin{array}{ccc}480& 0& 0\\ 0& 317& -6\\ 0& -6& 431\end{array}\right)\hspace{0.17em}\hspace{0.17em}\mathrm{N}\cdot \mathrm{m}$$

(77)

Both the values of matrices $Q_{11}^{}$, $Q_{22}^{}$ and $K_{\omega }$, $K_{\lambda }$ are rather large, which may lead to large control efforts for some domains of the initial conditions. However, the sets of these parameters are found in such a way that for any initial conditions in the domains (72) and (73), the control constraint (59) can not be violated. So, the control will always be less than feasible $u_{\max }=1\hspace{0.17em}\hspace{0.17em}\mathrm{N}\cdot \mathrm{m}$.

These two sets are used in the following section, where the algorithm operation is demonstrated.

The crucial proof for the stability is the fact that the control matrices and inertia tensor are diagonal in the same basis. The problem here is the fact that the inertia tensor used in the model can (and likely) differs from the one in the real system. Let $J_0$ be the nominal inertia tensor that is used in the mathematical model and $J=J_0\left(E_{3\times 3}+\epsilon j\right)$ be the real one, where $\epsilon$ is a small parameter that shows the difference between the nominal and real tensor and j is the symmetric matrix, the norm of which in some sense is one.

Consider matrix $K_{\omega }$ ($K_{\lambda }$ is analogous)

$$\begin{array}{l}{K}_{\omega }=R^{-1}{J}^{-1}{P}_{11}=R^{-1}{\left(E_{3\times 3}+\epsilon j\right)}^{-1}{J}_0^{-1}{P}_{11}\approx \\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\approx R^{-1}\left(E_{3\times 3}-\epsilon j\right)J_0^{-1}{P}_{11}=R^{-1}{J}_0^{-1}{P}_{11}-\epsilon R^{-1}j{J}_0^{-1}{P}_{11}.\end{array}$$

(78)

The first term is the positive definite matrix, while the second one in the general case can be even nonsymmetrical. However, the equivalent symmetric form matrix is

$${\omega }^T{R}^{-1}j{J}_0^{-1}{P}_{11}\omega =\frac{1}{2}{\omega }^T\left(R^{-1}j{J}_0^{-1}{P}_{11}+P_{11}^T{J}_0^{-1}j{R}^{-1}\right)\omega ,$$

(79)

so

$$K_{\omega }\approx R^{-1}{J}_0^{-1}{P}_{11}-\frac{1}{2}\epsilon \left(R^{-1}j{J}_0^{-1}{P}_{11}+P_{11}^T{J}_0^{-1}j{R}^{-1}\right).$$

(80)

The components of this matrix in the basis of principal axes of $J_0$ are

$$K_{\omega }=\left(\begin{array}{ccc}{k}_{11}^0-\epsilon k_{11}^1& -\epsilon k_{12}^1& -\epsilon k_{13}^1\\ -\epsilon k_{12}^1& k_{22}^0-\epsilon k_{22}^1& -\epsilon k_{23}^1\\ -\epsilon k_{12}^1& -\epsilon k_{23}^1& k_{33}^0-\epsilon k_{33}^1\end{array}\right).$$

(81)

Since $k_{ii}^0>0$, $k_{ij}^0$ and $k_{ij}^1$ have the same order while $\epsilon$ is a small parameter (the error is usually small), then the matrix $K_{\omega }$ is positive definite for the sufficiently small $\epsilon$.

7. Numerical Example

To demonstrate the typical system behavior, a numerical example is presented (Figure 3, Figure 4, Figure 5, Figure 6, Figure 7 and Figure 8). The system parameters for numerical simulation are presented in

Table 2. System parameters.

$m_s$,$\mathrm{kg}$	3500
$J_s$, $\mathrm{kg}\cdot {\mathrm{m}}^2$	$diag(\begin{array}{ccc}1.8& 0.9& 1.5\end{array})\times {10}^3$
$s_1$, $\mathrm{m}$	$\left(\begin{array}{ccc}0& 0& -0.75\end{array}\right)$
$s_2=-s_3$, $\mathrm{m}$	$\left(\begin{array}{ccc}0.5& 0& 0\end{array}\right)$
$m_1$, $\mathrm{kg}$	130
$J_1$, $\mathrm{kg}\cdot {\mathrm{m}}^2$	$\left(\begin{array}{ccc}8& 0& 0\\ 0& 7& -0.4\\ 0& -0.4& 2\end{array}\right)\times {10}^4$
$f_1$, $\mathrm{m}$	$\left(\begin{array}{ccc}0& 2.5& -12\end{array}\right)$
$m_{2,3}$, $\mathrm{kg}$	310
$J_{2,3}$, $\mathrm{kg}\cdot {\mathrm{m}}^2$	$\left(\begin{array}{ccc}180& -0.6& 0\\ -0.6& 730& 0\\ 0& 0& 910\end{array}\right)$
$f_2=-f_3$, $\mathrm{m}$	$\left(\begin{array}{ccc}4& 0& 0\end{array}\right)$

.

The numerical modelling is performed in the nonlinear model (10)–(14). Two sets of control parameters are taken. The simulation results are shown in the following figures.

Figures show that the control stabilizes the satellite and decreases the modal variable amplitudes. The process is rather slow since the control is small with respect to the total inertia tensor of the satellite with flexible elements. The peaks in Figure 5 and Figure 6 correspond to the switching between the sets of control parameters. As one can see from Figure 7 and Figure 8, this allows an increase in the convergence rate. The control level is almost three times less than the threshold. This is due to the fact that the approach guarantees (63) for each initial condition set in the domain (72). The evolution of the reaction wheel angular momentum in Figure 6 is due to the gravity gradient torque compensation.

8. Conclusions

The stabilization of a satellite with large flexible elements by means of reaction wheels only is shown in the paper. The stabilizing control based on the LQR is provided. The sufficient condition for the asymptotic stability is derived. This condition has an explicit form and can be checked once the spacecraft configuration is known. The choice of control parameters is based on the closed form solution of Riccati’s equation and the degree of stability of the system. The PSO usage allows one to find the parameters that give a rather good convergence rate and at the same time fulfill the control constraints.