Nonlinear Controller Design for a Flexible Spacecraft

Abstract

This paper combines the nonlinear Udwadia-Kalaba control approach with the Assumed Mode Method to model flexible structures and derives an attitude controller for a spacecraft. The study case of this paper is a satellite with four flexible cantilever beams attached to a rigid central hub. Two main topics are covered in this paper. The first one is the formulation of the equation of motion and the second one is the nonlinear controller design. The combination of these two techniques is able to provide a controller that damps the vibration of a flexible structure while achieving the desired rigid-motion state.

Introduction

Historically the flexibility of most satellites has not posed a major difficulty. This is because most of the flexible structures deployed so far have natural frequencies considerably higher than external disturbances and control actuation. Examples of cases where large structures have been modeled as rigid bodies can be found in [16] and [23]. In [15] and [3], even though rigid dynamics are used, the potential influence of the coupling between natural frequencies and external or internal disturbances is pointed out. This coupling is studied for two particular controllers over a large solar sail in [18]. Studies on how to model the deformation of large structures in space are also conducted in [4, 7, 22, 25].

In contrast to the historic perspective the interest of space industry in flexible structures has considerably increased during the last decades. See e.g. recent planetary society mission LightSail-2 [13] or DLR mission Tandem-L [11]. The current development in this area is driven by two main factors [24]: the miniaturization of satellites [1, 12, 14] and technological advances in related key areas, such as ultra-light deployable booms or sail films [2]. The main advantages of these structures are their low launching volume, as they are deployed once in orbit, and the low mass in relation to the final size of the structure. Examples of applications of this kind of technology are solar sailing [19], active deorbiting of debris [17] and solar power generation [9]. Such lightweight structures triggered the development of controllers with specific performance bounds [26] and fixed maneuver times [5]. Additionally, flexibility in combination with liquid propellant is becoming more and more relevant [6, 8].

This paper demonstrates a controller design methodology which takes into account the complex nonlinear motion of a flexible satellite. The chosen study case consists of a rigid body with four flexible cantilever beams attached to it, see Fig. 1.

Fig. 1

Study case: Satellite with four flexible cantilever beams

The objective is to develop a controller able to fulfill attitude commands and to damp the structural flexibility by applying forces and torques to the central body only. The paper is therefore organized as follows. Section “Equation of Motion Modeling” explains how the general equation of motion is formulated. The equations of motion are formulated by combining the Lagrange equation with the Assumed Modes Method (AMM). The AMM is utilized to convert the continuous variables, e.g. deformation of a beam, to a set of discrete variables, e.g. deformation of the beam in each eigenmode. Section “Model Implementation” particularizes the equations of Section “Equation of Motion Modeling” for the case of study.

Section “Control Approach” introduces the Udwadia-Kalaba [21] control approach, which uses the complete system model to compute the control command. It allows targeting with the controller all variables included in the description of the system, e.g. deformation of a beam in a particular eigenmode, and more complex variables that can be derived from those, e.g. potential energy due to deformation of several beams. Although the UK approach assumes full state knowledge and control, we demonstrate the possibility of reduced state control within the controller implementation in Section “Control Implementation”. Section “Controllers’ Performance” aims to study the performance of different controller designs. Finally, a conclusion sums up the results and provides some outlook to further studies.

Equation of Motion Modeling

The Udwadia-Kalaba control approach requires a model of the dynamics of the system. In this case, an analytical model of the equation of motion for the flexible satellites is built. The modeling approach used is based on the Lagrange equation in combination with the AMM, as described in [10, ch. 4.3.1]. This method is capable of generating meaningful models while maintaining a low order, which makes it suitable for control purposes.

Starting with the Lagrangian L = T − V, where T is the kinetic energy and V the potential energy, the generic equation of motion is expressed as shown in Eq. 1. Term q is the vector of generalized coordinates and Q represents all non conservative forces. Carrying out the derivation and distinguishing between the internal and external non conservative forces leads to the expression in the right where the internal non-conservative forces are expressed with the Rayleigh dissipation function D.

$$ \begin{array}{@{}rcl@{}} \frac{d}{dt}\left( \frac{\partial L}{\partial \dot{q}}\right)-\frac{\partial L}{\partial q}&=Q \Rightarrow \frac{d}{dt}\left( \frac{\partial T}{\partial \dot{q}}\right)-\frac{\partial T}{\partial q}+\frac{\partial V}{\partial q}+\frac{\partial D}{\partial \dot{q}}=Q_{ext} \end{array} $$

(1)

The expression for the kinetic energy (T), potential (V) and dissipation energy (D) can be expressed as

$$ \begin{array}{@{}rcl@{}} T =\frac{1}{2} \dot{q}^{T} M \dot{q} \qquad V = \frac{1}{2} q^{T} K q \qquad D = \frac{1}{2} \dot{q}^{T} F \dot{q}\qquad \end{array} $$

(2)

The term M states for the mass matrix, K for the spring constant matrix and F for the matrix of damping coefficients. These three matrices are square and have N × N terms, being N the length of the coordinates vector (q). These expressions allow performing directly the partial derivatives in Eq. 1 and arrive at the equation of motion (3), where Ω is the angular rate between inertial and body fixed reference frames. This expression fully describes the equation of motion and is used throughout the rest of this study.

$$ \begin{array}{@{}rcl@{}} \begin{array}{l} \overbrace{ M\ddot{q} + \dot{M}\dot{q} + {\varOmega} \times (M \dot{q}) }^{ \substack{ \frac{dx}{dt}_{A} = \frac{dx}{dt}_{B} + {\varOmega}_{BA}\times x \\ \Rightarrow \frac{d}{dt}\left( \frac{\partial T}{\partial \dot{q}} \right) = M\ddot{q}+\dot{M}\dot{q}+{\varOmega} \times (M \dot{q}) } } - \overbrace{ \frac{1}{2} \dot{q}^{T}\frac{\partial M}{\partial q}\dot{q} }^{\frac{\partial T}{\partial q}} + \overbrace{Kq}^{\frac{\partial V}{\partial q}} + \overbrace{F\dot{q}}^{\frac{\partial D}{\partial \dot{q}}} = Q \end{array} \end{array} $$

(3)

Having the equations of motion at hand the task is to find real valued expression for the matrices M, K and F. Therefore the Assumed Mode Methods (AMM) [10, ch. 4.3.1] is used. The AMM provides an approach to approximate the differential equation model of a continuous system and therefore its flexibility. To do this the deformation variable u is represented by a linear combination of time t dependent functions τ and space s dependent functions ϕ.

$$ u(s,t) = \sum\limits_{k=1}^{N} \tau_{k}(t)\phi_{k}(s) $$

(4)

The subindex k identifies the mode, being N the total number of modes considered. Each shape function ϕ_k describes a distribution of the deformation along the boom. The appropriate selection of these shape functions is one of the critical steps of this method. The terms τ_k contain the weights of each of these shape functions throughout time.

The considered study case accounts only for the bending (no torsion) of the (cantilever) beams that are part of the structure (see Fig. 1). The resulting number of coordinates for each boom is equal to the number of modes considered multiplied by the two orthogonal directions in which a beam can bend.

Combining Eqs. 3 with 4 allows integrating the flexibility into the equations of motion. The procedure to do this consists of the following steps

• 1. Describing the system and the system variables q to be taken into account.

• 2. Finding the expressions for the kinetic energy T, the potential V and the dissipation energy D as a function of those variables.

• 3. Redefining the continuous variables into sets of discrete variables using the AMM, assuming a suitable set of shape functions ϕ_k(s) can be found.

• 4. Finding matrices M, K and F.

• 5. Computing terms in Eq. 3.

Model Implementation

As previously mentioned, the study case consists of a satellite with four relatively long flexible cantilever beams, see Fig. 1. Here index I denotes the inertial reference frame and index c the body fixed reference frame (attached to the central part of the satellite). The mass properties, boom length etc. are given in Table 1.

Table 1 Structural properties

The state variables p to be studied are related to the rigid motion of the central hub of the satellite (r, v, 𝜃, ω) and to the flexible motion of each beam w.r.t. the central hub (u, $\dot {u}$).

$$ \begin{array}{@{}rcl@{}} p = \begin{bmatrix} r & \theta & u_{1} & u_{2} & u_{3} & u_{4} \end{bmatrix}^{T} \qquad \dot{p} = \begin{bmatrix} v & \omega & \dot{u}_{1} & \dot{u}_{2} & \dot{u}_{3} & \dot{u}_{4} \end{bmatrix}^{T} \end{array} $$

(5)

The next step is deriving expressions for the kinetic T, potential V and dissipation energy D as a function of p and $\dot {p}$. The total kinetic energy is the sum of the central part T_c and each beam T_bi (see Eq. 6). Here v_bi is the velocity at each point of the boom w.r.t. to the inertial reference frame, s is the longitudinal coordinate along the booms, index i denotes the boom and d the position of the point w.r.t. the body-fixed reference frame.

$$ \begin{array}{@{}rcl@{}} T = T_{c} + \sum\limits_{i=1}^{4} T_{bi} &=& \frac{1}{2} \omega^{T} I_{c} \omega + \frac{1}{2} m_{c} v^{2} + \sum\limits_{i=1}^{4}{{\int}_{0}^{L}} \frac{1}{2} \rho v_{bi}^{2} ds \\ v_{bi} &=& v + \omega \times d_{i} + \dot{u} \end{array} $$

(6)

Combining the previous expressions, the kinetic energy can be expressed as

$$ \begin{array}{@{}rcl@{}} T &=& \frac{1}{2} \omega^{T} I_{c} \omega + \frac{1}{2} m_{c} v^{T} v + 2 m_{b} v^{T} v + \frac{1}{2} \rho \sum\limits_{i=1}^{4} {{\int}_{0}^{L}} \left( - v^{T} [d_{i} \times] \omega \right.\\ &&\left.+ 2 v^{T} \dot{u} - 2 \dot{u}^{T} [d_{i} \times] \omega + \dot{u}^{T} \dot{u} + \omega^{T} [d_{i} \times]^{T}[d_{i} \times] \omega \right) ds \end{array} $$

(7)

With respect to potential and dissipation energy, only the terms related to the boom deformation are considered. The derivation of these equations (see expression (8)) can be found in [10].

$$ \begin{array}{@{}rcl@{}} V = \frac{1}{2} \sum\limits_{i=1}^{4} {{\int}_{0}^{L}} EI\left( \frac{ \partial^{2} u}{\partial s^{2}}\right)^{2}ds \qquad D = \frac{1}{2} \sum\limits_{i=1}^{4} {{\int}_{0}^{L}} k_{d}\left( \frac{ \partial^{2} \dot{u}}{\partial s^{2}}\right)^{2}ds \end{array} $$

(8)

It can be observed that the expressions derived for T, V and D contain continuous deformation variables u, which are dependent of time and space. In order to be able to define the matrices M, K and F (defined in Eq. 2), these variables need to be expressed as a sum of discrete coordinate. By utilizing the AMM space and time are decoupled and the integrals can be solved in the space domain beforehand.

$$ \begin{array}{@{}rcl@{}} u(s,t) = \sum\limits_{k} \tau_{k}(t)\phi_{k}(s) \rightarrow \int u(s,t) ds = \sum\limits_{k} \tau_{k}(t) \int \phi_{k}(s) ds \end{array} $$

(9)

As only bending deformations are considered and the section of the boom is assumed isotropic, only one set of shape functions ϕ_k is needed. These shape functions, taken from [10], describe the bending deformation of a cantilever beam and are expressed as

$$ \begin{array}{@{}rcl@{}} \phi_{k} (s) = 1-cos\left( \frac{k \pi s}{L}\right)+\frac{1}{2}(-1)^{k+1}\left( \frac{k \pi s}{L}\right)^{2} \end{array} $$

(10)

where k refers to the mode. The shape functions linked to the first five modes can be seen in Fig. 2. The decision to show the first five modes is purely due to visualization purposes. In this figure s is the normalized distance from the central part of the satellite to the boom tip. Utilizing the AMM, the state vector p can be redefined into q to become space independent.

$$ \begin{array}{ll} q = \begin{bmatrix} r & \theta & \tau_{1,x,1} &...& \tau_{i,j,k} &...& \tau_{4,z,N} \end{bmatrix}^{T} &\qquad i \in \{1,2,3,4\} \quad \text{beam index} \\ \dot{q} = \begin{bmatrix} v & \omega & \dot{\tau}_{1,x,1} &...& \dot{\tau}_{i,j,k} &...& \dot{\tau}_{4,z,N} \end{bmatrix}^{T} &\qquad j \in \{x,y,z\} \quad \text{coordinate} \\ &\qquad k \in \{1,...,N\} \quad \text{mode number} \end{array} $$

(11)

Based on this redefinition the final expressions for the matrices M, K and F can be assembled. Due to their considerable extension they are shown in Appendix I, in Eqs. 30 to 32.

Fig. 2

Shape functions ϕ_k(s) of the first 5 modes

The last open point is to find the remaining terms of the equation of motion (3). This includes $\dot {M}\dot {q}$, ${\varOmega } \times (M \dot {q})$, $\frac {1}{2} \dot {q}^{T}\frac {\partial M}{\partial q}\dot {q}$ and Q. The derivation of $\dot {M}$ is straightforward, as the only time-dependent variable appearing in M is related to the deformation of the beams in each mode.

The non-zero components of term ${\varOmega } \times (M \dot {q})$, which address the relative rotation between inertial and body-fixed reference frames, are:¹

$$ \begin{array}{@{}rcl@{}} \begin{array}{ll} {\varOmega} \times (M \dot{q})(1:3) &= \omega \times (M(1:3,1:3) v)\\ {\varOmega} \times (M \dot{q})(4:6) &= \omega \times (M(4:6,4:6) \omega) \end{array} \end{array} $$

(12)

The term $\frac {1}{2} \dot {q}^{T}\frac {\partial M}{\partial q}\dot {q}$ is more complex. In order to simplify its combination in the final expressions and avoid including tensors in the implementation, the term $\frac {1}{2} \dot {q}^{T}\frac {\partial M}{\partial q}$ is computed as a matrix. The resulting expressions are included in the Appendix I. The final step is to obtain the vector of generalized forces Q from the real forces that are applied to the system. This is done by calculating the virtual work that the real forces would perform as a function of the generalized coordinates q, as:

$$ \begin{array}{@{}rcl@{}} \delta W_{nc} = \sum\limits_{i} Q_{i} \delta q_{i} \end{array} $$

(13)

Equation 13 can be extended as δW = F_iδr_i + T₀δ𝜃₀, where F is force, T torque and 0 to 4 are the indices for the central hub and the four beam tips. The position is identified by r and the Euler angles by 𝜃. The control input is limited to F₀ and T₀. This finishes the definition of all terms in Eq. 3 such that a mathematical model for the flexible satellite in Fig. 1 is available. For further development the satellite system model itself was implemented in MATLAB^®; and validated against a classical FEM model.

Control Approach

We utilize the non-linear controller design approach proposed by Udwadia et al. in [21] and transfer it to our problem. The main principle of this control method is to formulate a constrained dynamic system, which follows the Gauss Principle of Least Constraint. The dynamic system itself consists of the plant dynamics, modeling constraints and all controller enforced constraints e.g. trajectory requirements. By applying the nowadays called Udwadia-Kalaba (UK) equation of motion [20, 27] the constraint force, which is equivalent to the control force, is directly visible and can be implemented within a classical control framework. In general this approach requires full state observability.

The following section summarizes briefly the non-linear control approach described in [21]. It starts with the generic equation of motion of a mechanical system. In the absence of constraints its equation of motion can be formulated as

$$ \begin{array}{@{}rcl@{}} M\ddot{q} = Q \Rightarrow \ddot{q} = M^{-1} Q = a \end{array} $$

(14)

Here M is a positive definite mass matrix, $\ddot {q}$ is the generalized acceleration a (the double time derivative of the state vector q) and Q the sum of all internal and external forces. Whenever there are arising additional constrains this generic equation of motion has to be rewritten by adding a constraint force component Q_c to Eq. 15.

$$ \begin{array}{@{}rcl@{}} M\ddot{q} = Q + Q_{c} \end{array} $$

(15)

The open problem is now to find a valid formulation for Q_c such that all constraints can be addressed and Eq. 14 holds. Such an expression is given in [21] by

$$ \begin{array}{@{}rcl@{}} Q_{c} = -K e \qquad K = M^{1/2} \left[ A M^{-1/2} \right]^{+} \qquad e = Aa - b \end{array} $$

(16)

where K can be interpreted as a feedback gain and term e as an error signal. Following this interpretation the force Q_c in Eq. 16 is the control force for the uncontrolled plant in Eq. 14.

From the control designer perspective we would like the plant in Eq. 14 to follow a specific guidance trajectory. This is interpreted within the UK framework as an additional constraint of one of the two forms

$$ \begin{array}{@{}rcl@{}} \phi_{i}(q,t) = 0 \qquad \psi_{i}(q,\dot{q},t) = 0 \end{array} $$

(17)

When a guidance is existing as holonomic (ϕ) or as nonholonomic (ψ) expression it can be differentiated once or twice w.r.t. time to obtain the relation

$$ \begin{array}{@{}rcl@{}} A(q,\dot{q},t)\ddot{q} = b(q,\dot{q},t) \end{array} $$

(18)

which gives the remaining terms A and b in the control force of Eq. 16.

In a real world system it is often hard to directly apply trajectory constraints of the form Eq. 17. They require that the state variable q is already contained in the described manifold. Therefore, it is more suitable to generalize the constraints to

$$ \begin{array}{@{}rcl@{}} \ddot{\phi} + {\varSigma} \dot{\phi} + {\varGamma} \phi = 0 \qquad \dot{\psi} + {\varLambda} \psi = 0 \end{array} $$

(19)

which ensures that over time the trajectory can converge but does not require that it matches exactly right from the beginning. The variables Σ, Γ and Λ are constant tuning parameters, which describe the desired convergence rates.

When this equation of motion formulation is followed cost function J is automatically minimized at each instant of time.

$$ \begin{array}{@{}rcl@{}} J(t) = Q_{c}(q,\dot{q},t)^{T} \ M^{-1}(q,t) \ Q_{c}(q,\dot{q},t) \end{array} $$

(20)

This cost function is representing the constraint (control) force which is selected by nature for any constrained motion. It is also possible to generalize the cost function, which gives some freedom for the control designer.

$$ \begin{array}{@{}rcl@{}} J(t) = Q_{c}(q,\dot{q},t)^{T} \ N(q,t) \ Q_{c}(q,\dot{q},t) \end{array} $$

(21)

The use of this different “mass” matrix N means that the resulting controller is not “equivalent” to Nature’s. The matrix N needs to be chosen carefully, as it has a great impact in the performance of the controller and can even limit the coordinates over which an actuation can be exerted. According to [21] the closed form solution of the control gain in Eq. 16 changes with the usage of this new weighting matrix N to.

$$ \begin{array}{@{}rcl@{}} K = N^{-1/2} A^{+}_{N} \qquad A_{N} = A M^{-1} N^{-1/2} \end{array} $$

(22)

A summary of the different ways in which the parallelism between constrained and controlled systems can be noticed is presented in Table 2.

Table 2 Equivalences of constrained and controlled systems [21]

Control Implementation

The general nonlinear control law is given by the control force Q_c, which is a combination of Eqs. 16 and 22.

$$ \begin{array}{@{}rcl@{}} Q_{c} = -N^{-1/2} \left[AM^{-1} N^{-1/2} \right]^{+} (Aa - b) \end{array} $$

(23)

For our problem Q_c represents the forces and torques that are acting on the structure of the satellite. The following paragraphs identify all terms inside (23) and explain what they are used for.

Unconstrained System $\leftrightarrow$ Uncontrolled Plant - a

Here we extract from the general equation of motion (3) the unconstrained equation of motion (14) to identify the acceleration a. This acceleration term is defined by the physics of the problem and can not be adjusted by the control designer.

$$ \begin{array}{@{}rcl@{}} a = \ddot{q} = M^{+}Q = M^{+} \left( - \dot{M}\dot{q} - R + \frac{1}{2} \dot{q}^{T} \frac{\partial M}{\partial q} \dot{q} - Kq - F\dot{q} \right) \end{array} $$

(24)

Constraints $\leftrightarrow$ Trajectory Requirements - A, b

The first desired control target is to drive the angular rate towards a constant guidance value ω_g. This is a nonholonomic constraint of the same form as Eq. 17, where ω is equal to $\dot {q}$. To express this constraint within the UK framework we need its first derivative and make use of Eq. 19, which allows a convergence toward the desired angular rate. The constraint equation is therefore:

$$ \begin{array}{@{}rcl@{}} \psi_{\omega}(q, \dot{q}, t)=0 \quad \Rightarrow\quad \omega -\omega_{g} = 0 \quad \Rightarrow\quad \dot{\omega}+{\varLambda}_{\omega}(\omega-\omega_{g}) = 0. \end{array} $$

(25)

Another constraint targets the time derivative of the boom deformations $\dot {u}$. The idea is to bring all motion within the booms down to zero, which is a holonomic constraint. Based on the AMM, this constraint can be rewritten for each mode and boom. The generalized coordinate is $\dot {q}$ equals $\dot {\tau }$. By taking its 1st derivative w.r.t. time in order to allow for some convergence we end up with expression (26). Such a formulation allows formulating requirements for each single boom and each single mode. It is even possible to control only a individual mode motion of a single boom.

$$ \begin{array}{@{}rcl@{}} \phi_{\tau}(\dot{q}, t) = 0 \quad \Rightarrow\quad \dot{u} = 0 \quad \Rightarrow\quad \dot{\tau} \phi = 0 \quad \Rightarrow\quad \ddot{\tau} + {\varSigma}_{\tau} \dot{\tau} = 0 \end{array} $$

(26)

The third constrain formulation addresses a general vibration damping by minimizing the dissipation energy. This formulation has basically the same effect as the previous constrain but much simpler to evaluate at run time. To use this formulation we need to derive it w.r.t. time once again.

$$ \begin{array}{@{}rcl@{}} \psi_{F}(q, \dot{q}, t)= 0 \quad \Rightarrow\quad \dot{q}^{T} F \dot{q} = 0 \quad \Rightarrow\quad 2\dot{q}^{T}F\ddot{q}+ {\varLambda}_{F} \dot{q}^{T}F\dot{q} = 0 \end{array} $$

(27)

Now we have all constraint (control goal) equations in place to assemble the matrix A and vector b of Eq. 18.

$$ \begin{array}{@{}rcl@{}} A = \begin{bmatrix} 0_{[3,3]} & I_{[3,3]} & 0_{[3,n]} \\ 0_{[n,3]} & 0_{[n,3]} & I_{[n,n]} \\ 0_{[1,3]} & 0_{[1,3]} & 2\dot{q}^{T} F \end{bmatrix} \qquad b = \begin{bmatrix} - {\varLambda}_{\omega} (\omega - \omega_{g}) \\ - {\varSigma}_{\tau} \dot{\tau} \\ - {\varLambda}_{F} \dot{q}^{T} F \dot{q} \end{bmatrix} \end{array} $$

(28)

It should be noted here that the last two constrains are nearly equivalent and it does not make sense to implement them both at the same time. For the case under study we used different combinations of the matrices A, and b, where we skipped the rows depending on the irrelevant trajectory requirements.

Controller Tuning - N

The “modified mass” matrix N represents the cost, see Eq. 21, of exerting an actuation over each of the generalized coordinates. It allows to tune the desired controller behavior.

$$ \begin{array}{@{}rcl@{}} N &= \begin{bmatrix} m_{T} I_{[3,3]} & 0 & 0 \\ 0 & I^{-1} & 0 \\ 0 & 0 & 10^{15} I_{[n,n]} \end{bmatrix} \end{array} $$

(29)

The weighting terms related to the forces at the central satellite hub are set equivalent to the total mass m_T. The terms related to the corresponding torques are set to the inverse moment of inertia I^− 1. Additional control actuation should not take place. Therefore, the terms of coordinates linked to the beam deflection are penalized very strong and set to 10¹⁵. This could potentially lead to a badly conditioned matrix N, which could give some problems during its inversion. Therefore, this number was some kind of trial and error to figure out acceptable numerical stability and proper control force calculation.

Controllers’ Performance

This section demonstrates the potential use of the previously defined controllers over the structure defined in Table 1. In these initial results, for the sake of simplicity, only the first flexible mode is considered. Additionally, the decision was made to limit the control force to 0.01 N and torque to 0.01 Nm (in each axis, exerted in the central part of the structure). No extensive tuning of the controllers was conducted, setting Λ_F and Λ_ω and Σ_τ to 1. The definition of the external disturbances, torque and (distributed) force, can be found in Appendix II.

Four different test scenarios are shown in Table 3. The 1^st test is conducted in the absence of a controller and servers as basic reference. The 2^nd contains a controller that focuses on damping the vibration of the booms, the 3^rd on controlling the angular rate and the 4^th combines both control objectives. The remaining of this section contains analyses of this four test scenarios in both time and frequency domains.

Table 3 Controllers implemented

Time Domain

Two analyses are conducted, focusing on the response to a force/torque impulse and to an external disturbances profile. In relation to the response to impulses, two study cases are considered: 1) response to an application of a torque of 1Nm for 1s in the Z-axis and 2) response to an application of a force of 1N for 1s in the Z-axis. The results, evolution of total energy, are presented in Fig. 3.

Fig. 3

Total energy responses to torque and force impulse in the Z-axis

For the response to a torque impulse, in the absence of control (1) the energy remains constant after the impulse (t = 5s). Controller 2 reduces the energy by damping the vibration of the booms but does not reduce the angular rate, leading to a constant energy. Cases 3 and 4 show the same behavior, both reducing indefinitely the energy. This means that controller 3 also damps the vibration of the booms, which is only true when the vibrations are visible in the angular motion.

The response to a force impulse does not result in an angular motion. Therefore, controller 2, which aims to damp the angular rate, shows the same performance as no controller (case 0). In this case the energy of the system is related to the velocity of the structure and the deformations of the booms. Both controller 2 and 4 are able to damp the vibrations of the booms. The velocity is not among the control variables, so there is some remaining energy.

After analyzing the response to force/torque impulses, the performance of the different controllers is studied for a force/torque profile of external disturbances. The profile used can be found in Appendix II. The initial state is given in Table 4.

Table 4 Case of study setup

The evolution of the system is studied for the 4 different cases defined in Table 3. In order to give an overall idea of the state of the system, the total energy is not enough. Therefore, the evolution of angular rate, potential energy and dissipation energy is shown. Figure 4 shows the response for scenarios 1 and 2. It can be observed that in scenario 1 both V, D and ω increase with time.

Fig. 4

Evolution of V, D and ω (X-axis blue, Y red, Z yellow), Scenarios 1 & 2

In scenario 2 the vibration is damped considerably fast after each change in the external disturbances. The order of magnitude of V and D is at least one order of magnitude lower than without a controller. The angular rate is not controlled because the controller is only targeting the boom vibrations.

Scenarios 3 and 4 are shown in Fig. 5. In scenario 3, it can be seen that the controller is able to reach and maintain the desired ω. In the evolution of V and D, it can be seen that the vibrations in the booms are similar to those appearing in the absence of control (scenario 1). However, the oscillations in ω are damped as the angular rate error reduces. Before that reduction the control command surpasses the maximum actuation capability and the actuator is saturated.

Fig. 5

Evolution of V, D and ω (X-axis blue, Y red, Z yellow), Scenarios 3 & 4

In scenario 4 two periods can be distinguished, before and after the ω error is low enough for the actuators not to be saturated. In the first period the control command is cut and most of the damping of the structure is lost. In the second period the controller is able to both maintain the desired ω and considerably damp the system.

Frequency Domain

The initial state vector is defined in Table 4. The term ω_g is set to zero here. Due to the non-linear character of the system, this analysis was conducted as follows. First an external disturbance with a particular frequency is defined (e.g. torque around the z-axis at 1 Hz). Then the simulation is conducted until the response of the system is periodic. Finally, the amplitude of the signals for the study variables is measured. An exemplary output of such an input is shown in Fig. 6 for the angular rate.

Fig. 6

Evolution of ω_y under a torque in y-axis with a frequency of 1.9744 rad/s, Scenario 1

The inputs considered are forces and torques in X,Y and Z-axis, applied in the central part of the satellite. The outputs studied are V, D and ω. The range of frequencies studied is between 1 and 30 rad/s, which includes the natural frequencies of the system. For better readability Fig. 6 in Appendix III shows all numerical determined input-output transfer functions over the frequency grid. The natural frequencies of the system (frequencies where there is a peak in response) are clearly observable, both w.r.t. external forces and external torques.

The response of ω to a torque in the same axis is similar for the different controllers implemented, particularly around the natural frequencies of the system (≈ 12 rad/s). The controllers are able to reduce the peak response almost one order of magnitude. At lower frequencies control controller 2 has a lower performance, as it only targets the angular acceleration indirectly as a consequence of damping the system. Considering the response of V and D to a input torque, the three controllers show again a similar performance, reducing the peak response by a factor of 40. At low frequencies, the effect of the controller in the potential energy is lower, while the effect in the dissipation energy is maintained. This is because it is the dissipation energy the variable targeted by the controllers.

When the input considered is a force, the effect of controller 3 is negligible. The other control approaches (2 and 4) have similar performances, leading to a peak reduction of factor 50 in both the potential and dissipation energy. Similarly, to the case with torques as inputs, at low frequencies the effect of the controller in the potential is lower, while the effect in the dissipation energy is maintained. It is also noticeable that the natural frequency of the system is shifted to the left, from 12 to 2.6 rad/s.

Conclusions

This study shows the potential of using the Udwadia-Kalaba control approach in combination with a method to model flexible structures. This control approach has two main advantages: 1) it allows implementing multiple complex guidance instructions and 2) maintains the non-linear dynamics of the system. The method used to model the structure is based in the Lagrange equation combined with the AMM (assumed modes method), and provides a trade-off between accuracy and simplicity of the model. It integrates the flexibility of the structure into its equation of motion, which can be used directly to derive the controller.

Among the main limitations that are foreseen w.r.t. the use of this control approach in this kind of structures are: 1) the computational cost of each control step, 2) the high frequency at which the controller could have to run and 3) the need of estimating the deformation of the booms. From a modeling perspective, the main difficulty lays in finding shape functions able to accurately represent the flexible modes of the structure.

The results show that controllers derived using this approach are able to actively damp the structure, reducing the dissipation energy and stabilizing the potential of the structure, while controlling the angular rate.

The directions of future research in this area include, but are not restricted to, the following:

• Further tuning of the controllers, including constraints and cost function.

• Simplified experimental tests of the control approach.

• Detailed model verification.

• Study of additional constraints, such as limitations in the torque direction.

• Study of the potential applicability to other structures.

Appendix

I. Equation of Motion

This appendix contains the complete expression of the terms in the equation of motion (3).

The mass matrix is described by Eq. 30. Here subindex i refers to the axis, j to the boom and k to the mode. The term m is the mass of the central part of the satellite (CS) and m_b of one boom (B), I_u is an unitary 3x3 matrix, ρ the linear density of B, L the length of B, [dx] is a matrix containing the cross product by the position of a point of B (see expressions Eqs. 6 and 7), ds is a differential in the longitudinal axis of B, Φ_{i, j} is shown in expression (34), and I is the moment of inertia of CS. The spring constant and damping coefficient matrices are described in Eqs. 31 and 32, respectively. Here EI is the bending stiffness of B and k_d the damping coefficient.

The expressions related to the shape functions ϕ_k(s) are shown in Eq. 35, where δ_kl is a Kronecker delta. These expressions are used for particularizing the equations of motion for the shape functions chosen.

$$ M = \begin{bmatrix} (m+4m_{b})I_{u}&{\rho{\int}_{0}^{L}}\sum [dx]_{j}ds&{\rho{\int}_{0}^{L}} {\Phi}_{1,1}&{\rho{\int}_{0}^{L}}{\Phi}_{1,2}&...&{\rho{\int}_{0}^{L}}{\Phi}_{4,N-1}&{\rho{\int}_{0}^{L}}{\Phi}_{4,N}\\ -{\rho{\int}_{0}^{L}}\sum [dx]_{j}ds &I+{\rho{\int}_{0}^{L}}-\sum [dx]_{j}^{2} ds&{\rho{\int}_{0}^{L}} [dx]_{1}^{T}{\Phi}_{1,1}ds&{\rho{\int}_{0}^{L}} [dx]_{1}^{T}{\Phi}_{1,2}ds&...&{\rho{\int}_{0}^{L}} [dx]_{4}^{T}{\Phi}_{4,N-1}ds&{\rho{\int}_{0}^{L}} [dx]_{4}^{T}{\Phi}_{4,N}ds\\ {\rho{\int}_{0}^{L}} {\Phi}_{1,1}&{\rho{\int}_{0}^{L}} {\Phi}_{1,1}[dx]_{1}ds&{\rho{\int}_{0}^{L}} {\Phi}_{1,1}^{2}& {\rho{\int}_{0}^{L}} {\Phi}_{1,1}{\Phi}_{1,2} & ...& 0 & 0\\ {\rho{\int}_{0}^{L}} {\Phi}_{1,2}&{\rho{\int}_{0}^{L}} {\Phi}_{1,2}[dx]_{1}ds&{\rho{\int}_{0}^{L}} {\Phi}_{1,1}{\Phi}_{1,2}& {\rho{\int}_{0}^{L}} {\Phi}_{1,2}^{2} & ...& 0 & 0\\ ...&...&...& ... & ... & ... & ...\\ {\rho{\int}_{0}^{L}} {\Phi}_{4,N-1}&{\rho{\int}_{0}^{L}} {\Phi}_{4,N-1}[dx]_{4}ds&0& 0 & ...& {\rho{\int}_{0}^{L}} {\Phi}_{4,N-1}^{2} & {\rho{\int}_{0}^{L}} {\Phi}_{4,N-1}{\Phi}_{4,N}\\ {\rho{\int}_{0}^{L}} {\Phi}_{4,N}&{\rho{\int}_{0}^{L}} {\Phi}_{4,N}[dx]_{4}ds&0& 0 & ...& {\rho{\int}_{0}^{L}} {\Phi}_{4,N-1}{\Phi}_{4,N}& {\rho{\int}_{0}^{L}} {\Phi}_{4,N}^{2} \\ \end{bmatrix} $$

(30)

$$ K = \begin{bmatrix} 0&0&0&0&...&0&0\\ 0&0&0&0&...&0&0\\ 0&0&{{\int}_{0}^{L}} EI(\frac{ \partial^{2} {\Phi}_{1,1}}{\partial s^{2}})^{2}ds& {{\int}_{0}^{L}} EI\frac{ \partial^{2} {\Phi}_{1,1}}{\partial s^{2}}\frac{ \partial^{2} {\Phi}_{1,2}}{\partial s^{2}}ds & ...& 0 & 0\\ 0&0&{{\int}_{0}^{L}} EI\frac{ \partial^{2} {\Phi}_{1,1}}{\partial s^{2}}\frac{ \partial^{2} {\Phi}_{1,2}}{\partial s^{2}}ds&{{\int}_{0}^{L}} EI(\frac{ \partial^{2} {\Phi}_{1,2}}{\partial s^{2}})^{2}ds & ...& 0 & 0\\ ...&...&...& ... & ... & ... & ...\\ 0&0&0& 0 & ...& {{\int}_{0}^{L}} EI(\frac{ \partial^{2} {\Phi}_{4,N-1}}{\partial s^{2}})^{2}ds & {{\int}_{0}^{L}} EI\frac{ \partial^{2} {\Phi}_{4,N-1}}{\partial s^{2}}\frac{ \partial^{2} {\Phi}_{4,N}}{\partial s^{2}}ds\\ 0&0&0& 0 & ...& {{\int}_{0}^{L}} EI\frac{ \partial^{2} {\Phi}_{4,N-1}}{\partial s^{2}}\frac{ \partial^{2} {\Phi}_{4,N}}{\partial s^{2}}ds&{{\int}_{0}^{L}} EI(\frac{ \partial^{2} {\Phi}_{4,N}}{\partial s^{2}})^{2}ds\\ \end{bmatrix} $$

(31)

$$ F = \begin{bmatrix} 0&0&0&0&...&0&0\\ 0&0&0&0&...&0&0\\ 0&0&{{\int}_{0}^{L}} k_{d}(\frac{ \partial^{2} {\Phi}_{1,1}}{\partial s^{2}})^{2}ds& {{\int}_{0}^{L}} k_{d}\frac{ \partial^{2} {\Phi}_{1,1}}{\partial s^{2}}\frac{ \partial^{2} {\Phi}_{1,2}}{\partial s^{2}}ds & ...& 0 & 0\\ 0&0&{{\int}_{0}^{L}} k_{d}\frac{ \partial^{2} {\Phi}_{1,1}}{\partial s^{2}}\frac{ \partial^{2} {\Phi}_{1,2}}{\partial s^{2}}ds&{{\int}_{0}^{L}} k_{d}(\frac{ \partial^{2} {\Phi}_{1,2}}{\partial s^{2}})^{2}ds & ...& 0 & 0\\ ...&...&...& ... & ... & ... & ...\\ 0&0&0& 0 & ...& {{\int}_{0}^{L}} k_{d}(\frac{ \partial^{2} {\Phi}_{4,N-1}}{\partial s^{2}})^{2}ds & {{\int}_{0}^{L}} k_{d}\frac{ \partial^{2} {\Phi}_{4,N-1}}{\partial s^{2}}\frac{ \partial^{2} {\Phi}_{4,N}}{\partial s^{2}}ds\\ 0&0&0& 0 & ...& {{\int}_{0}^{L}} k_{d}\frac{ \partial^{2} {\Phi}_{4,N-1}}{\partial s^{2}}\frac{ \partial^{2} {\Phi}_{4,N}}{\partial s^{2}}ds&{{\int}_{0}^{L}} k_{d}(\frac{ \partial^{2} {\Phi}_{4,N}}{\partial s^{2}})^{2}ds\\ \end{bmatrix} $$

(32)

$$ \begin{array}{@{}rcl@{}} \dot{q}^{T}(t)\frac{\partial M(q,t)}{\partial q}(a,:) \!&=&\! \frac{1}{2} \!\!\begin{bmatrix} ({\rho{\int}_{0}^{L}}\frac{\partial\sum [dx]_{j}}{\partial q_{ijk}} ds) \omega\\ (-{\rho{\int}_{0}^{L}}\frac{\partial\sum [dx]_{j}}{\partial q_{ijk}} ds) v - ({\rho{\int}_{0}^{L}} \frac{\partial\sum [dx]_{j}^{2} }{\partial q_{ijk}}ds) \omega + {\sum}_{lmn}(({\rho{\int}_{0}^{L}}\frac{\partial [dx]_{l}^{T}}{\partial q_{ijk}} {\Phi}_{l,m}ds)\dot{q}_{nml})\\ ({\rho{\int}_{0}^{L}} {\Phi}_{1,1}\frac{\partial [dx]_{1}}{\partial q_{ijk}}ds)\omega\\ ({\rho{\int}_{0}^{L}} {\Phi}_{1,2}\frac{\partial [dx]_{1}}{\partial q_{ijk}}ds)\omega\\ ...\\ ({\rho{\int}_{0}^{L}} {\Phi}_{4,N}\frac{\partial [dx]_{4}}{\partial q_{ijk}}ds)\omega\\ \end{bmatrix}^{T}\\ a &=& 6+i+3(k-1)+3N(j-1);\quad i=x,y,z; k=1...N; j=1...4; \end{array} $$

(33)

$$ {\Phi}_{j,k} = \begin{bmatrix} \phi_{xjk}& 0&0\\ 0&\phi_{yjk}&0\\ 0&0&\phi_{zjk} \end{bmatrix} $$

(34)

$$ \begin{array}{@{}rcl@{}} {{\int}_{0}^{L}} \phi_{k} ds &=& L \left( 1+\frac{1}{6}(-1)^{k+1} (k \pi)^{2}\right); \quad {{\int}_{0}^{L}} s \phi_{k} ds \\ &=& L^{2}\left( \frac{1}{2}+ \frac{1}{8}(-1)^{k+1}(k \pi)^{2}-((-1)^{k}-1)\frac{1}{(k\pi)^{2}}\right) \end{array} $$

(35)

$$ \begin{array}{@{}rcl@{}} {{\int}_{0}^{L}} \phi_{k} \phi_{l} ds &=& L\left( 1+(-1)^{k+l}\left( \left( \frac{k}{l}\right)^{2}+\left( \frac{l}{k}\right)^{2}\right)+\frac{1}{2}\delta_{kl} +(-1)^{k+1} \frac{k^{2} \pi^{2} }{6}\right.\\ &&\left.+(-1)^{l+1}\frac{l^{2} \pi^{2} }{6}+(-1)^{k+l}\frac{k^{2}l^{2}\pi^{4}}{20}\right) \end{array} $$

(36)

$$ \begin{array}{@{}rcl@{}} {{\int}_{0}^{L}} \phi^{\prime\prime}_{k} \phi^{\prime\prime}_{l} ds &=& \left( (-1)^{k+l}+\frac{1}{2}\delta_{kl}\right)\frac{k^{2}l^{2}\pi^{4}}{L^{3}}; \quad \phi_{k} (L) = 1+(-1)^{k+l}\\ &&+\frac{1}{2}(-1)^{k+l}(k^{2}\pi^{2}); \quad \phi^{\prime}_{k} (L) = (-1)^{k+l}\frac{k^{2}\pi^{2}}{L} \end{array} $$

(37)

II. External Disturbances

Fig. 7

External disturbances distribution

Fig. 8

Profile of external disturbances

Fig. 9

Evolution of control actuation (X-axis blue, Y red, Z yellow), Scenario 4

The external disturbances consist in a distributed force and a torque applied in the center of the satellite. The force is distributed between the center of the satellite (50%) and the tip of each boom (12.5%) (see Fig. 7). This distribution is chosen to model a external pressure such as atmospheric drag or solar pressure.

The evolution of the external disturbances for the study case is defined as shown in Fig. 8. Here the left diagram shows the torques and the right one the forces applied to the satellite plant.

Figure 9 shows the control actuation applied in scenario 4, Fig. 5, Section 1. It can be seen that the force command focuses only on damping the vibration while the torque command aims also to achieve a certain angular rate ω.

III. Response in the Frequency Domain

Fig. 10

Response of the different controllers in the frequency domain