ADCS Design for a Sounding Rocket with Thrust Vectoring

This paper addresses the development of an attitude determination and control system (ADCS) for a sounding rocket using thrust vector control (TVC). To design the ADCS, a non-linear 6 degrees-of-freedom (DoF) model for the rocket dynamics and kinematics is deduced and implemented in simulation environment. An optimal attitude controller is designed using the linear quadratic regulator (LQR) with an additional integral action (LQI), and relying on the derived linear, time-varying, state-space representation of the rocket. The controller is tested in the simulation environment, demonstrating satisfactory attitude tracking performance, and robustness to model uncertainties. A navigation system is designed, based on measurements available on-board, to provide accurate real-time estimates on the rocket’s state and on the aerodynamic forces and moments acting on the vehicle. These aerodynamic estimates are used by an adaptive version of the controller that computes the gains in real time after correcting the state-space model. Finally, the ADCS is the result of the integration of the attitude control and navigation systems, with the complete system being implemented and tested in simulation, and demonstrating satisfactory performance.


Introduction
The main motivation behind this work is the development of an Attitude Determination and Control System (ADCS) for a future sounding rocket from a student rocketry team from Instituto Superior Técnico (IST), named Rocket Experiment Division (RED). The ADCS design assumes that the rocket uses Thrust Vector Control (TVC) technology as the actuation method, and aims to control the rocket's pitch and yaw angles. The roll angle is assumed to be controlled by an additional roll control system whose design is out of the scope of this work.
During the atmospheric flight phase of a rocket, stabilisation can be achieved through the use of aerodynamic fins. With a correct design of the fins, the vehicle can be made naturally stable [1]. However, the rocket is subjected to various external disturbances, such as wind gusts, which prevent the vehicle to follow a desirable, pre-calculated trajectory or, even more intense, completely destabilise it [2]. It is then clear the necessity of having an active attitude control and stabilisation system that not only ensures the stability of the rocket, but allows to actively correct its trajectory in order to achieve specific mission goals. As for the actuation method, Thrust Vector Control (TVC), or thrust vectoring for short, is used by most launch vehicles and works by redirecting the thrust vector in order to create a control torque [3]. With respect to other actuation techniques, like actively controlled fins, TVC allows for a wider range of operating conditions and provides better efficiency [4].
The control system design tends to be very conservative in the aerospace industry. Restricting the dynamic analysis to accommodate more sophisticated control design techniques risks the later realisation that such restrictions would have to be lifted and would invalidate the control design. Among the classical techniques, the Proportional-Integral-Derivative (PID) control is on the core of most commonly used launch vehicle control systems [3,5]. Although widely used, PID control has its downsides when it comes to robustness and external disturbances rejection [6]. The problem of controlling ascending launch vehicles is dominated by parameter uncertainty, which in face of the lack of robustness of the 1 3 PID controller may be a concerning issue. Moreover, the rocket flight parameters considerably change throughout the flight. To overcome this, gain scheduling techniques have been proposed, that rely on the linearization of the dynamics at different operating conditions. Still in the linear domain, the use of optimal controllers, such as the Linear Quadratic Regulator (LQR), provides more robustness and ensures a locally optimal solution for a given cost function [6,7]. As a way to improve the robustness of linear time-varying controllers, real-time parameter estimators can be introduced in the control loop to form an adaptive control system. The online identification of system parameters allows the controller to act on a more accurate representation of the system dynamics [8].
Non-linear control and estimation techniques have also been proposed [9,10], and come with the advantage of ensuring a global solution, not dependent on the specific mission nor vehicle. However, this type of control and estimation laws often have to simplify complete non-linear dynamics in order to obtain a global solution. If relevant dynamics are discarded, the system might fail in a real implementation scenario. Besides that, these methods all have particular design characteristics which make it harder to develop a standardised verification and validation procedure to meet the imposed system requirements [3,5].
Although several solutions to the rocket attitude control problem can be found in the literature [7,11,12], many fail to capture all the relevant dynamics and/or oversimplify the problem, while most assume full-state knowledge, creating a considerable gap between theoretical design and implementation. Hence, the main contribution of this paper is a robust ADCS, which integrates both the navigation and control systems, relying on computationally efficient algorithms and that can be implemented in sounding rockets through readily available low-cost components. Furthermore, the design process considers the 6 DoF, as opposed to restricting the analysis to the pitch plane, as well as the time-varying nature of the system, focussing on the entire trajectory rather than a single operating point, contrarily to what is found in the literature when using linear optimal control and estimation techniques. In this way, the implementation of the system in a real case scenario is facilitated.
To achieve this goal, several intermediate contributions were necessary, always having in mind the reliability and robustness of the proposed solutions: an original generic state-space representation for linear and optimal control design in the 6 DoF; a gain-scheduled optimal pitch and yaw controller resorting to the LQR technique with added integrative action; and a navigation system, based on measurements available on-board, to provide accurate estimates on the rocket's state, including a novel linear time-varying parameter estimator to estimate in real-time aerodynamic forces and moments acting on the vehicle.

Rocket Dynamics and Kinematics Modelling
To design the ADCS, a mathematical model that represents the translational and rotational dynamics and kinematics of the rocket in the 6 DoF is necessary.

Assumptions
Some assumptions are used to derive the model. The rocket is considered to be a rigid body, meaning no elastic behaviours are modelled. This assumption is considered valid for control system design given the smaller size of typical sounding rockets, and consequent reduced impact of elastic behaviour on the overall dynamics. The rocket is assumed to be axially symmetric, as well as the mass allocation, which means that the principal inertia axes coincide with the body axes, the centre of mass is on the longitudinal axis, and the aerodynamic behaviour is identical in both the pitch and yaw planes. This assumption is standard given that launch vehicles, and more specifically sounding rockets, are designed in order to respect this symmetry. Finally, neither the curvature nor rotation of the Earth are taken into account, which is also a reasonable assumption considering the typical altitude and ground distance covered by the class of vehicles under study.

Reference Frames
To describe the dynamics and kinematics of the rocket, it is crucial to define the reference frames to be used. Two reference frames are used: a body-fixed one (Fig. 1a), where the equations of motion are written; and an inertial space-fixed one (Fig. 1b). The body-fixed reference frame has its origin located in the centre of mass of the vehicle. The x-axis ( X b ) is along the rocket's longitudinal axis, while the z-axis ( Z b ) and y-axis ( Y b ) complete the orthogonal reference frame. As for the inertial space-fixed reference frame, given that neither the curvature nor the Fig. 1 Body-fixed a and inertial b reference frames rotational motion of the Earth are taken into account, a simple orthogonal frame centred in the launch location is used. The x-axis ( X e ) is pointing upwards, so that for a zero inclination launch the x-axes of both reference frames are aligned; and the other two axes ( Y e and Z e ) are preferably aligned with a pair of cardinal directions.
With the reference frames detailed, it is necessary to define the coordinate transformation between them. This is done using a sequential rotation of the body frame relative to the Earth frame defined by the three Euler angles where is the Euler angle of rotation of the body around the x-axis of the Earth frame, also known as roll; is the Euler angle of rotation of the body around the y-axis of the Earth frame, also known as pitch; and is the Euler angle of rotation of the body around the z-axis of the Earth frame, also known as yaw. The Euler angles describe the attitude of the rocket, representing the variables to be controlled by the attitude control system. The coordinate transformation from the body frame to the Earth frame is then defined by the following transformation matrix [13]: where c and s stand as abbreviations for the trigonometric functions. The inverse transform, from the Earth frame to the body frame, is defined by the transpose ( R T ).

External Forces and Moments
From the dynamic point of view, sounding rockets experience four main forces during a flight: Weight, Thrust, and the aerodynamic forces-lift and drag.

Gravity Model
Considering the Earth as a perfect sphere, the gravitational acceleration, g, is assumed to vary only with altitude. This variation is given by where g 0 is the gravitational acceleration constant at surface level, R E is the mean Earth radius, and h is the altitude. The gravitational force, expressed in the Earth frame, is given by E F g = −mg 0 0 T , where m is the rocket's mass.
To obtain it in the body frame, the rotation matrix R T is used,

Propulsion Model
The propulsion model was derived using equations mainly obtained from [4], considering ideal propulsion and all its underlying assumptions. The thrust produced by the rocket motor is simply where ṁ is the mass flow rate, v e is the effective exhaust velocity, p e is the nozzle exit pressure, p a is the atmospheric pressure, and A e is the nozzle exit area. Two separate contributions can be identified: the dynamic one, caused by the exhaust of the expanded combustion gases; and the static, caused by the pressure gradient between the nozzle exit and the atmosphere. Considering that the most common propulsion technology for sounding rockets is solid propulsion, and that it is the technology used by RED, a model that uses the internal combustion equations, available in [4,14], as well as the solid propellant characteristics, is implemented to calculate the thrust produced by the motor, T, the associated mass flow rate, ṁ , and the nozzle exit pressure, p e . The propellant considered for this work was the mixture of potassium nitrate with sorbitol (KNSB), a propellant commonly used in student rocketry known as "rocket candy", with its properties present in [15]. The mass and geometry of the propellant grains can be altered to obtain different thrust profiles.

TVC Actuation
By controlling the direction of the thrust force (or vector), TVC actuation produces torques that act on the rocket's centre of mass, influencing its rotation in pitch and yaw. The decomposition of the propulsive force in the three body axes can be done as illustrated in Fig. 2.
According to it, the thrust vector is decomposed using the angles p and y , where p is the gimbal angle that, on its own, produces a pitching moment, and y is the one that produces a yawing moment. Using these angles, the propulsive force in the body frame is given by B F p = T cos p cos y − T cos p sin y − T sin p T , a n d t h e c o r r e s p o n d e n t c o n t r o l m o m e n t i s B M = 0 − T sin p l T cos p sin y l T [6], where l is the moment arm, which corresponds to the distance between the nozzle gimbal point and the centre of mass of the rocket.

Aerodynamic Forces and Moments
The rocket will be subjected to aerodynamic forces and moments resulting from its interaction with the fluid medium composing the atmosphere. Starting by the forces, they are expressed in the body axes according to where C A is the axial aerodynamic force coefficient, C Y is the lateral aerodynamic force coefficient, C N is the normal aerodynamic force coefficient, q is the dynamic pressure and S is a reference area, usually corresponding to the cross sectional area of the fuselage. The axial and normal aerodynamic forces correspond to the body axes components of lift and drag, and are related through the aerodynamic angles-the angle of attack, = arctan w rel ∕u rel , and the sideslip angle, = arcsin v rel ∕V rel , where u rel , v rel , and w rel are the components of the relative velocity vector with respect to the atmosphere, and V rel its magnitude. The force coefficients can be determined using a linear relation with the aerodynamic angles, C Y = C Y and C N = C N , whose derivatives ( C Y and C N ) depend mainly on the angle itself and Mach number.
As for the aerodynamic moments, in the body axes they are given by C l , C m , and C n are, respectively, the rolling, pitching, and yawing moment coefficients, and d is a reference length, usually corresponding to the diameter of the fuselage. If the reference moment station is defined as the centre of pressure, and its location, x cp , can be determined, the reference moments are zero and the moment coefficients take the form where the static stability margin, SM = (x cp − x cm )∕d , intuitively appears, p, q, and r, are the body angular velocities, and C l p , C m q , C ṁ , C n r , and C ṅ are all aerodynamic damping coefficients.

Translational Motion
By applying Newton's second law, and taking into account that the body frame is a rotating one, we obtain the translation dynamics, where v = [u v w] T is the velocity vector in the body frame, = p q r T is the angular velocity vector in the body frame, S(.) is a skew-symmetric matrix, and the mass derivative term has been included in the propulsive force. By substituting the external forces in (1), the dynamics can be particularised in the body acceleration components:

Rotational Motion
Euler's equation for rigid body rotational motion yields where J is the inertia matrix. Following the axial symmetry assumption, the cross-products of inertia can be assumed as zero and the y and z terms can be assumed equal, resulting in a diagonal matrix, J = diag (J l , J t , J t ) , where J l denotes the longitudinal inertia and J t denotes the transverse inertia. By substituting the inertia matrix J and the external moments in where r represents the rolling moment caused by the additional roll control system and r accounts for external disturbances. Given that the aim of this work is to control the pitch and yaw angles only, the additional roll control system is assumed to be able to reject disturbances on this axis ( r ) and control the roll angle. Its inclusion in the model is only for the sake of completeness and is not considered during the design. Finally, the rotational kinematics are given by the time derivative of the Euler angles [13]: It is noted that using the Euler angles a singularity arises for = ± 2 , however, the way the reference frames are defined prevents the rocket to reach this attitude inside the admissible range of operation (far from horizontal orientation). By grouping (2), (4), and (5) the 6 DoF non-linear model of the rocket is fully defined.

Reference Rocket
A specific rocket is needed to serve as reference for the ADCS design. In this way, a preliminary design for a future RED's rocket with Thrust Vector Control is performed. The rocket is designed to have a burning phase coinciding with the full duration of the climb, so that TVC can be used to control its attitude up to apogee. It is also required that the terminal velocity is inside a safe range to allow the correct activation of the recovery system. To meet these design goals, the solid motor parameters are iteratively tested using the propulsion model, and the flight for a vertical undisturbed trajectory is simulated. Tables 1 and 2 respectively present the main rocket characteristics and the simulation results.

Model Linearization
To design an optimal, linear controller, it is necessary to obtain a linear version of the model and respective statespace representation. The non-linear model can be linearized at equilibrium points of the system using a Taylor series expansion, considering small perturbations. For the case of a rocket, conditions change considerably throughout the flight, hence, it is not correct to choose a single equilibrium point to linearize the system. Instead, a reference trajectory is selected and the system is linearized at multiple operating points. The outcome is a linear time-varying system. When obtaining the linear version of the system, it is advantageous to consider some assumptions: the roll rate, p, is considered to be zero as it will be controlled by an external roll control system, reducing the order of the system by one; the wind is assumed to be zero, allowing to directly use the linear velocity in the body frame in the aerodynamic terms; the actuator dynamics are not included in the model; and the system parameters are considered to be constant at the linearization points, removing the dependencies on the state variables when computing the Taylor derivatives. By applying a Taylor series expansion to the non-linear system around the operating points, and considering these assumptions, a linear time-varying system in the perturbation domain is obtained, that can be represented in the state-space form: where A(t) and B(t) are the state-space matrices given by the first-order Taylor derivatives with respect to system states and inputs, respectively, calculated at the operating points.
Regarding the attitude reference that defines the reference trajectory, a varying pitch trajectory, in which the controller restricts the motion to the pitch plane (yaw equal to zero) and makes the rocket deviate from the vertical to later recover it, is selected. In this way, it is ensured that the apogee is reached further away from the launch site, increasing safety. Figure 3 shows the reference pitch angle and rate over time.
(6a)  Since the system is naturally unstable, it is necessary to find the time evolution of the nominal control inputs that allows the rocket to nominally follow the trajectory, defined by an attitude reference over time. This is done using a PID controller that in simulation, without perturbations, is able to stabilise the vehicle and track the attitude reference. The input values over time are then stored to use as predetermined feedforward control inputs.
It is possible to identify two distinct sections of the reference trajectory: a first section up to t = 25 s in which motion is strictly vertical, and a second section up to burnout in which pitch is varying. For the varying pitch section, we have that: 0 = 0 = 0 , v 0 = 0 , r 0 = 0 , and y 0 = 0 . This results in a simplified version of the statespace representation, for which the longitudinal and lateral modes are decoupled: In the vertical section, we have that: 0 = 0 = 0 = 0 , v 0 = w 0 = 0 , q 0 = r 0 = 0 , and p 0 = y 0 = 0 . This results in a simplified version of the decoupled state-space representation, for which u and are no longer states of the system. Finally, it is important to determine the location of the system poles throughout the nominal trajectory to derive the open-loop stability. For a time-varying system, the stability is not mathematically guaranteed with this method, however, the study is carried out to understand the behaviour of the system throughout the flight. Figure 4 details the pole evolution (from blue to green) during the vertical section and the poles at t = 50 s, which serves as example for the distribution type during the varying pitch section. Looking at the different pole distributions during the flight, some conclusions can be made. First, the system has poles located in the right-hand side of the complex plane for the entire trajectory, meaning that it is naturally unstable. This was expected due to the negative static stability margin caused by the absence of aerodynamic fins. Second, during the vertical section, there is an equivalence between the lateral and longitudinal modes, verified by the identical pole distributions, which is due to the symmetry of the vehicle and to the equality in the nominal values of the corresponding states. Thirdly, the system has complex conjugate poles with positive real part during the first seconds of the flight, which indicates natural unstable oscillatory behaviour, after which all poles start to be located on the real axis. Finally, it is concluded that the velocity of the rocket is a driving factor for the response of the system-as velocity increases during the flight, the system is seen to have higher magnitude poles and hence faster response. This is attributed to the fact that at higher velocities the magnitude of the aerodynamic forces and moments is also higher, causing higher accelerations on the system when the inputs are actuated.

Linear Quadratic Integral (LQI) Control
Using the linear time-varying state-space representation of the system, a linear quadratic regulator (LQR) is designed with the addition of an integral action, also known as linear quadratic integral control (LQI). The LQR is a technique that finds the optimal gain matrix k for the linear control law u = −k x , which minimises a quadratic cost function given by where Q is a positive semi-definite matrix and R is a positive definite matrix [16]. In the cost function, the quadratic form ′ Q represents a penalty on the deviation of the state x from the origin, and the term ′ R represents the cost of control, making Q and R the tuning parameters for the resultant controller. Using the infinite-horizon version, which means taking T as infinity, the solution which minimises the cost function and guarantees closed-loop asymptotic stability is the constant gain matrix k = R −1 B T P , where P is the solution to the Algebraic Riccati Equation (ARE) P A + A T P − P B R −1 B T P + Q = 0.
Since the system is time-varying, the ARE has to be solved for models coming from each linearization point, resulting in a set of gain matrices to be selected, or scheduled, throughout the flight.
The LQR feedback control law ideally drives the states of the system in the perturbation domain to zero, ensuring that the nominal values throughout the trajectory are followed. However, it does not guarantee a zero tracking error for non-zero references in terms of attitude. In order to have a zero reference tracking error, and to increase the robustness of the controller, an integral action that acts on the attitude tracking error is added, according to the scheme in Fig. 5.
Let the difference between the reference signal, r , and the output of the system, y , (the tracking error) be the time derivative of the state-space variables that result from adding the referred integrator, . The state-space representation of the resulting regulator can be obtained by combining the open-loop state-space representation with the feedback law, where z = x T is the augmented state vector and C is the output matrix that selects the output of the system from the original state vector ( y = C x ). The optimal gain K is obtained by solving the ARE using the rearranged system matrices, Considering the decoupling between the longitudinal and lateral modes, the decoupled augmented state vectors are z lon = [ u w q i ] T and z lat = [ v r i ] T , where i and i are the integral states. This implies that the A , B are divided into the longitudinal and lateral modes, and that the C matrix for the lateral mode is the one that selects the yaw angle, while for the longitudinal mode is the one that selects the pitch angle.
The design degree of freedom is the selection of the tuning matrices Q and R , which will also be divided into First of all, setting all non-diagonal entries to zero, and only focussing on the diagonal ones, allows for a more intuitive matrix selection given by the "penalty" method. According to this method, the diagonal entries of the Q matrix will determine the relative importance of the state variables in terms of origin tracking performance, while the diagonal entries of the R matrix allow to directly adjust the control effort for each input. Therefore, the weighting matrices have the following generic format, separated for each mode, Given the nature of the TVC actuation, trying to control the linear velocities would conflict with the attitude control, specially for non-zero attitude references. Hence, the linear velocity related terms are set as zero. By doing this, the associated gains will have negligible magnitude, allowing to use partial state feedback with k lon = k q k k i and k lat = k r k k i . The tuning parameters are iteratively adjusted looking at the closed-loop poles and at the step response performance in the linear domain, including the actuator dynamics, modelled as a first-order system. Regarding the closed-loop poles, the control law allowed to stabilise all operating points, placing all closed-loop poles in the lefthand side of the complex plane. Table 3 details the step response parameters for multiple operating points.

Navigation System Design
So far, it was assumed that the control system has access to an exact full-state measurement. In reality, it is necessary to have a navigation system, composed by sensors and estimators, capable of providing an accurate estimate on the state vector. For the case of rockets, and taking into account the state variables to measured, it is common to use an Inertial Measurement Unit (IMU), composed by accelerometers, gyroscopes, barometers, and magnetometers, and a Global Navigation Satellite System (GNSS) receiver. The estimator architecture was based on [17]. It is composed by three main filters and a pre-processing unit (PU), according to the scheme in Fig. 6.
The pre-processing unit (PU) combines the magnetometer and accelerometer readings, m r and a r , to obtain an indirect measurement on the Euler angles, r . Then, the first filter is an Attitude Complementary Filter (ACF), which will use the Euler angles readings and the measured angular rates from the gyroscopes, r , to provide a filtered attitude estimate, ̂ , and an estimate on the bias of the three angular rates, b , to correct the signal from the sensor. The second one is a Position Complementary Filter (PCF), which merges the position readings from the GNSS receiver, translated into the inertial frame, p r , and the acceleration measurements from the accelerometer to provide an estimate on the velocity vector, v . This filter is also self-calibrated since it accounts for the bias in the three acceleration readings, b a . Finally, a Linear Parameter Estimator (LPE), uses the control inputs, velocities, angular rates and attitude pre-filtered values to give a final estimate on the state vector, x , and parameters, ̂ .

ACF
For the ACF, it is assumed that the Euler angles measurement is corrupted by Gaussian white noise, w , as well as the angular rates readings, w , and that the gyroscope bias is described by a constant term with additional Gaussian white noise. Considering this, the filter is based on the kinematic equations for the Euler angles (5), using directly the Euler angles readings in the process matrices to allow for the use

PCF
For the PCF, both the position and acceleration measurements are considered to be corrupted by Gaussian white noise, w p and w a , and the accelerometer bias is also described by a constant term with additional Gaussian white noise. This filter is also kinematic, considering the following equations of motion: where p is the position in the inertial frame and a is the acceleration expressed in the body frame. The state-space representation of the filter is then obtained, where stands for the identity matrix and for the matrix of zeros, both of dimension 3 by 3. The rotation matrix R is calculated using the Euler angles estimate from the ACF ( ̂ ). The individual gain matrices , and , each with dimension 3 by 3, can once again be computed considering the vertical attitude time-invariant to define the rotation matrix R , so as to obtain time-invariant Kalman gains.

LPE
The robustness of the LQR is limited since the controller is designed considering a nominal evolution of model parameters that might considerably differ from the real evolution during the mission. Amongst the model parameters, the ones related with the aerodynamic properties of the rocket are subjected to an higher level of uncertainty, due to the difficulty in obtaining accurate aerodynamic coefficients and derivatives of the rocket for a broad range of velocities and aerodynamic angles. In this way, an online parameter estimator is proposed so that the controller acts on an informed value of the aerodynamic parameters. The aerodynamic parameters are hidden under the aerodynamic force and moment coefficients. Since a first estimate on these quantities is available using the stored aerodynamic data, a proportional error factor is multiplied in each aerodynamic force and moment, corresponding to the parameters to be estimated.
The estimator design follows along the methodology proposed in [18], where an hovercraft control system is designed based on dynamic parameters identification, which details a generic parameter estimator for time-varying systems, linear in the parameters.
The previously detailed rocket model (section 2.4) is rearranged by including the proportional error factors, a x , a y , a z , m , and n on the aerodynamic forces and moments, making B = −q C A S a x q C Y S a y −q C N S a z T and B = 0 q C m S d m q C n S d n T , where the aerodynamic rolling moment is discarded due to the additional roll control system. After substituting the rearranged aerodynamic forces and moments in the rocket model, and considering the linearity in the parameters to be estimated, the nonl i n e a r d i f fe r e n t i a l e qu a t i o n s t a ke t h e fo r m ̇x = f(x, t) + G(x, t) , where x = u v w q r T and = a x a y a z m n T . Using state augmentation with the parameter vector, , and assuming full-state measurements, y , are available, this system can be written in statespace form as in which the full-state measurement assumption allows to regard the system as linear, and the parameters are assumed to be slowly varying. The G(y, t) and f(y, t) matrices are easily obtained using the derived rocket model with the inclusion of the correction factors, and are not here presented to improve readability. In order to design the estimator for this system, it is necessary for it to be observable. In the reference, it is demonstrated that the system is observable if and only if there exists no unit vector d , with the dimension of the parameter vector, such that ∫ t t 0 G(y, ) d ⋅ d = 0 . Taking the time derivative in both sides and substituting for the rocket dynamics, the equivalent non-observability condition is where d i , for i = 1, 2, 3, 4, 5 , are the components of the unit vector, and the simplification is due to m, J, q , d, and S being always different from zero. It is possible to infer that the system is observable only when the aerodynamic force and moment coefficients are all different from zero, since if one of them is not, the unit vector with d i = 1, where i corresponds to component multiplying the null coefficient, satisfies the non-observability condition. However, given a null coefficient, the correspondent correction factor is the unobservable parameter, meaning that estimates can still be obtained for the remaining ones. Nevertheless, to ensure full observability, if the pre-calculation of a given coefficient results in zero, it can be forced to a small non-zero value. After verifying that the system can be made observable, a Kalman filter represents a simple and easily tunable solution for the estimation of the system state, resulting in the following state-space representation for the LPE:

Adaptive LQI Control
Resorting to the real-time estimates on the aerodynamic error coefficients, the LQI controller gains can be computed on-board instead of scheduling the pre-calculated gains for each operating point. This is done by rewriting the statespace representation with the inclusion of the estimated parameters, and solving the ARE on-board with the updated state-space models. The rearranged system dynamics matrix, A , can be easily achieved and is not here presented.

Simulation Environment
To test and validate the proposed ADCS in the complete non-linear model, a realistic simulation model is implemented in MATLAB/Simulink ® environment. The model is composed by several subsystems in order to completely transcribe the derived dynamics and kinematics, generate the environmental properties, and compute the model-varying parameters.
The environmental properties are generated by the atmospheric, wind, and gravitational models. The 1976 U.S. standard atmosphere was implemented, which describes the evolution of temperature and pressure with altitude using average annual values, from which density and speed of sound are derived. Wind is introduced through the summation of the average horizontal wind components from the U.S. Naval Research Laboratory horizontal wind model with wind gusts added from the Dryden model, both available as Simulink blocks. Finally, the gravitational model is implemented according to the equations in Sect. 2.3.1.
Several varying model parameters have to be computed during simulation. The ideal thrust force and mass flow rate are predetermined using the propulsion model detailed in Sect. 2.3.2, and the static, atmospheric pressure-dependant thrust component is added during the simulation. The aerodynamic properties, i.e. the aerodynamic derivatives and centre of pressure location, are stored in look-up tables and are selected according to the instant values of the aerodynamic angles and Mach number. The mass properties are also computed during the simulation, including the mass, inertia, and centre of mass, which vary due to the propellant consumption.
The equations of motion used in the simulation environment are the ones presented in Sect. 2.4. It is important to note that some assumptions were used when deriving the model and, although considered valid for design, can have impact on the expected performance, obtained in simulation, when in a real case scenario. Elastic modes might be excited by the control action if the associated frequencies are similar, causing undesired oscillatory behaviour; asymmetries may cause the centre of mass to be dislocated from the x-axis of the body, which imposes additional effort on the control action; and non-linear aerodynamic effects may cause unexpected behaviour, as well as unaccounted effects caused by the rotation of the Earth, such as the Coriolis acceleration.

ADCS Parameters
In this section, some of the parameters used in the simulation environment are detailed. These include the model used for the actuators, the control system gains, and the covariance matrices obtained for the navigation system.

Actuators' Model
The actuators' dynamics are modelled using a continuous time first-order transfer function for each input ( p and y ), considering a servo-actuated system. The transfer function is given by where r is the actuator angular response and is the time constant. In addition, servo motors normally have a saturation value for the rotation velocity, which can be modelled by a rate limiter block in Simulink. The time constant and angular velocity limit values were retrieved from typical high grade servo motors, and are equal to 0.02 s and 1 full rotation per second, respectively.

Control Gains
The tuning of the Q and R matrices for each mode yielded the time evolution for the controller gains throughout the nominal trajectory detailed in Fig. 7.
The gains remain approximately constant given that the tuning matrices were left constant for all operating points, except for the ones associated with the longitudinal mode during the varying pitch section, which were tuned in order to reduce the control effort and avoid saturation.

Estimation Covariances
Each of the individual filters composing the navigation system follows the Kalman filter structure [19], meaning that the tuning parameters are their respective covariance matrices, Q and R . The Q corresponds to the covariance of the process noise, and R corresponds to the covariance of the measurement noise. The R matrices can be derived by referring to the noise properties of the on-board sensors, while the Q matrices are iteratively adjusted looking at the simulation results in the realistic environment, yielding r = 1 s + 1 Q acf = 10 −1 , Q pcf = 10 −2 , Q lpe = 10 −2 10 , R acf = , R pcf = 4 , R lpe = 10 .

Navigation System
The navigation system is tested and is able to reject the noise introduced by the sensors, remove the bias, and provide an accurate estimate on the state of the rocket. Figure 8 presents the pitch angle and pitch rate estimation by the ACF, while Fig. 9 presents the crossrange position and longitudinal velocity estimation by the PCF, both as exemplification of the performance of the system. The linear parameter estimator is also tested in simulation, by inducing errors in the aerodynamic coefficients, and is able to correctly estimate the parameters. Figure 10 presents the results for a simulation in which errors were induced in the aerodynamic coefficients associated with pitch plane motion ( C A , C N , and C m ). The affected parameters ( a x , a z and m ) are seen to correctly converge to the expected values and the associated estimation errors for the aerodynamic forces and moments are minimised.

LQI Control
The LQI controller is implemented in the simulation model and tested by adding wind, with and without gusts, as external perturbation. Table 4 displays the results in terms of attitude tracking performance and control effort, not only for the LQI controller, but also for a tested PID controller for comparison.
It is noted that the LQI controller provides better attitude tracking for the same control effort with respect to the PID. In the yaw plane, the results for the PID are significantly worse since it is very affected by the initial wind perturbation. The step response is also analysed (Tab. 5).
Once again the LQI displays satisfactory performance, close to the design values (Tab. 3), and significantly better than the classical PID. Additionally, a robustness analysis is carried out, in which the model parameters are varied in percentage. For the assumed parameter uncertainties, the LQI controller shows high robustness. (c)

Adaptive LQI Control
Due to the high robustness of the non-adaptive LQI controller, the adaptive version is not able to produce significant performance improvements.

Complete ADCS
The complete ADCS is tested by integrating the attitude control and navigation systems. Table 6 details the attitude tracking performance and the control effort with wind gusts present, in comparison with the results for the control system alone without sensor noise.
As expected, there is a performance decrease. However, it is still satisfactory.

Sensitivity Analysis
Finally, a sensitivity analysis was performed to determine the robustness of the system to model uncertainties. Several system parameters, including dry mass, inertia, Thrust, centre of mass position, and aerodynamic coefficients, were altered independently, inside admissible ranges in terms of percentage of the original value: The system showed sufficient robustness, being able to stabilise the plant for all the variations under study. The parameter which demonstrated the highest influence in the performance of the control system was the position of the centre of mass ( x cm ), with the results shown in Table 7.
A lower value, meaning a position closer to the tip of the rocket, causes the moment arm for the thrust vector actuation to be higher, which increases the control authority. At the same time, the natural instability of the rocket reduces. In this way, the tracking performance increases when the centre of mass moves closer to the tip, while the control effort decreases.

Conclusions
With the conclusion of this work, it is possible to state that the primary goal has been achieved: the successful design of an attitude determination and control system applicable to sounding rockets with thrust vectoring. The design process was described in a generic way to ensure that the system can be easily applied to different vehicles under the same category. Nevertheless, the future implementation of the system in a student-built sounding rocket was always taken into account, as it was the initial motivation behind this work.
As future work, it would be of interest to develop nonlinear controllers for the attitude control problem in order to compare the performance of said controllers with the developed ones. Particularly, the designed linear parameter estimator could be used for a non-linear control system that  requires accurate information on the aerodynamic forces and moments to guarantee its correct functioning. Moreover, both the developed simulation model and navigation system can be verified and validated using real flight data from sounding rockets launched by RED. In fact, the simulation environment shall be improved by including phenomena yet to be modelled, such as elastic modes, non-linear aerodynamic effects, the curvature and rotation of the Earth, and body asymmetries, to further verify the system before implementation in a real case scenario. Finally, RED is currently developing small-scale prototypes to test the TVC technology and the associated navigation and control systems. In this way, it is intended to implement the techniques in here developed to such prototypes and to analyse all the results coming from test campaigns.

Conflict of interest
The authors declare that they have no known conflict of interest.
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:// creat iveco mmons. org/ licen ses/ by/4. 0/.