Nonlinear Controller Design for a Flexible Spacecraft

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 ultralight 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.
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" Fig. 1 Study case: Satellite with four flexible cantilever beams 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.
The expression for the kinetic energy (T), potential (V) and dissipation energy (D) can be expressed as 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.
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 φ.
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. 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,u).
The next step is deriving expressions for the kinetic T , potential V and dissipation energy D as a function of p andṗ. 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.
Combining the previous expressions, the kinetic energy can be expressed as 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].
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.
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 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  The non-zero components of term Ω × (Mq), which address the relative rotation between inertial and body-fixed reference frames, are: 1 The term 1 2q T ∂M ∂qq is more complex. In order to simplify its combination in the final expressions and avoid including tensors in the implementation, the term 1 2q T ∂M ∂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: Equation 13 can be extended as δW = F i δr i + T 0 δθ 0 , 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 0 and T 0 .
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 Here M is a positive definite mass matrix,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.
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 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 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 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 tö 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.
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.
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.
A summary of the different ways in which the parallelism between constrained and controlled systems can be noticed is presented in Table 2.

Control Implementation
The general nonlinear control law is given by the control force Q c , which is a combination of Eqs. 16 and 22.
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 ↔ 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.
Constraints ↔ 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 toq. 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: Another constraint targets the time derivative of the boom deformationsu. 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 isq equalsτ . 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.
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.
Now we have all constraint (control goal) equations in place to assemble the matrix A and vector b of Eq. 18.
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.
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 15 . 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.

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   Fig. 3 Total energy responses to torque and force impulse in the Z-axis 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.
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.
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.  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.
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. 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.
-Study of additional constraints, such as limitations in the torque direction.
-Study of the potential applicability to other structures.
Acknowledgements Open Access funding enabled and organized by Projekt DEAL. We appreciate the constructive comments and suggestions made by the reviewers.

Compliance with Ethical Standards
Conflict of interests The authors declare that there is no conflict of interest regarding the publication of this paper.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommonshorg/licenses/by/4. 0/.  , 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.

Generals
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.      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 ω.