DEVELOPMENT OF AN ACTIVE CONTROL

SYSTEM FOR A CANARD CONTROLLED

SOUNDING ROCKET

Theodore Youds

School of Mechanical Engineering

University of Leeds

May 2023
 MECH5080M Team Project

Author(s) Name: Theodore Youds

SID: 201230040

Supervisor(s) Name: Dr Jongrae Kim

Examiner Name: Dr Greg de Boer

Date of Submission: 03-05-2023

Development of an Active Control System

for a Canard Controlled Rocket

ii
SCHOOL OF MECHANICAL ENGINEERING

MECH5840
MECH5080M TEAM PROJECT 45 credits

TITLE OF PROJECT

Development of an Active Control System for a Canard Controlled Rocket

PRESENTED BY Theodore Youds

OBJECTIVES OF PROJECT

1. Create a mathematical model of the system dynamics from first

principles.
2. Implement this model in Simulink in order to predict the behaviour of
the system.
3. Design a controller that can deflect the canards in order to maintain
more vertical flight.

IF THE PROJECT IS INDUSTRIALLY LINKED TICK THIS BOX

×
AND PROVIDE DETAILS BELOW

COMPANY NAME AND ADDRESS:

Collins Aerospace

Wolverhampton

WV10 7EH

THIS PROJECT REPORT PRESENTS OUR OWN WORK AND DOES NOT CONTAIN ANY
UNACKNOWLEDGED WORK FROM ANY OTHER SOURCES.

INDUSTRIAL MENTOR:
SIGNED Theodore Youds DATE 3rd May 2023

iii
 Table of Contents

Acknowledgements ...................................................................................................vi

Abstract .................................................................................................................... vii

Nomenclature .......................................................................................................... viii

Symbols ............................................................................................................... viii

Acronyms ............................................................................................................... x

1 Introduction......................................................................................................... 1

1.1 Project Aptos ............................................................................................... 1

1.2 Individual Aims and Objectives .................................................................... 1

1.3 Report Structure .......................................................................................... 1

2 Literature Review and Theoretical Background .................................................. 2

2.1 Sounding Rocket Stability ............................................................................ 2

2.2 Closed Loop Control Systems ..................................................................... 3

2.3 Robustness Analysis ................................................................................... 3

3 Mathematical Modelling ...................................................................................... 4

3.1 Frames of Reference ................................................................................... 4

3.2 Pitch and Yaw Moments .............................................................................. 4

3.2.1 Pitch / Yaw Corrective Moment ............................................................ 4

3.2.2 Pitch / Yaw Damping Moment .............................................................. 5

3.2.3 Pitch / Yaw Canard Moment ................................................................. 5

3.3 Roll Moments .............................................................................................. 6

3.3.1 Roll Damping Moment .......................................................................... 6

3.3.2 Roll Canard Moment............................................................................. 7

3.4 Equation of Motion Formulation ................................................................... 7

3.5 Nonlinear Modelling ..................................................................................... 8

4 Closed Loop Controller Design ........................................................................... 9

4.1 States, Inputs, and Disturbances ................................................................. 9

4.2 State Space Formulation ............................................................................. 9

4.3 LQR Design ............................................................................................... 10

iv
 4.4 Gain Scheduling ........................................................................................ 11

4.5 Servomotor Position Control ...................................................................... 12

4.6 Robustness Analysis ................................................................................. 13

4.7 Controller Performance ............................................................................. 14

5 Conclusion ....................................................................................................... 15

5.1 Achievements and Impact ......................................................................... 15

5.2 Future Work............................................................................................... 15

6 References ....................................................................................................... 16

Appendix A – Supporting Formulae Derivation ........................................................ 19

Normal Force Coefficient ..................................................................................... 19

Calculating 𝑪𝑵𝜶 for Axisymmetric Body Components .......................................... 20

Calculating 𝑪𝑵𝜶 for Fins ...................................................................................... 20

CP for Body Components .................................................................................... 22

CP for Fins ........................................................................................................... 22

Full Rocket CP ..................................................................................................... 22

Appendix B – Linearisation ...................................................................................... 23

Roll ...................................................................................................................... 23

Pitch..................................................................................................................... 23

Yaw...................................................................................................................... 24

State Space Formulation ...................................................................................... 25

Appendix C – Expected Parameter Variations ......................................................... 26

Appendix D – Canard Moment Diagrams................................................................. 27

Appendix E – Enhanced Wind Modelling Results .................................................... 28

v
 Acknowledgements

Firstly, I would like to acknowledge my Aptos teammates Ben, Joe, Liam, and Jordan.
The achievements of this project speak for themselves, and it has been a pleasure
working with you. Next, I would like to thank our supervisor Dr Jongrae Kim for his
support of the project, as well the other academic staff within the School of Mechanical
Engineering, especially Dr Greg de Boer and Dr Bilal Kaddouh.

We wouldn’t have been able to conduct any flight testing without approval and
mentoring from the UK Rocketry Association. Thank you to Colin Rowe, Chris Brown,
and Paul Carter for your advice and guidance with all of LURA’s rocketry endeavours.
I am also extremely grateful for John Bonsor up at the Fairlie Moor Rocketry Site in
Scotland. John has always been an advocate for LURA and allowed us to use his
range to conduct preliminary flight tests in April 2023. His friendly, caring nature helped
make both launch days an incredibly enjoyable process.

Finally, I would like to thank my girlfriend Katy, my parents, my housemates, friends,

and family. You have supported me through everything over the past four years at the
University of Leeds and I wouldn’t have been able to do this without you. Thank you.

vi
 Abstract

Project Aptos aimed to develop a modular, canard controlled active stability system for
a sounding rocket. This report focuses specifically on the development of the control
algorithms used to command the system. An introduction to the project aims and
objectives is given, specifically relating to the key stakeholders in the project: the Leeds
University Rocketry Association (LURA) and Collins Aerospace. Relevant literature is
discussed and existing work on control systems for orbital rockets, sounding rockets,
and missile systems is drawn upon.

A three degree of freedom rotational model of a canard-controlled sounding rocket is

derived and subsequently modelled using MATLAB and Simulink. These equations of
motion are linearised using first order Taylor series techniques and small angle
approximations, and a linear quadratic regulator (LQR) system is designed to control
the rocket’s attitude during ascent.

A separate PID controller is implemented to improve the performance of servomotor

position control such that it meets the system requirements. The MATLAB System
Identification app is then used to fit a transfer function model to the servomotor step
response data, allowing for accurate modelling within Simulink.

The control system’s performance is analysed, and it is proven that performance is

sufficient for test flights to commence. The system’s robustness is investigated over a
range of expected variations in the rocket’s geometry using a Monte Carlo approach.
Finally, the Chernoff bound is applied in order to prove 95% accuracy and confidence
in the system’s performance and stability robustness.

vii
 Nomenclature

Symbols

𝛼 Angle of attack
𝐴𝑟𝑒𝑓 Reference area (rocket cross-sectional area)
𝐵 Body frame of reference
𝛽 Sideslip angle
𝐶1 Corrective moment coefficient
𝐶2 Damping moment coefficient
𝐶𝑑𝑝 Roll damping moment coefficient derivative w.r.t roll rate
𝐶𝑁𝛼0 Normal force coefficient derivative of fin 2D aerofoil
𝐶𝑁𝛼 Normal force coefficient derivative w.r.t. angle of attack
𝐶𝐿 Lift coefficient
𝐶𝐿𝛼 Lift coefficient derivative w.r.t angle of attack
𝑑𝑏𝑜𝑑𝑦 Body tube diameter
𝐝 Disturbance vector, {𝛼 𝛽} 𝑇
𝛿 Robustness analysis confidence = 1 − 𝛿
𝜖 Robustness analysis accuracy
𝐹 Canard lift force, e.g. 𝐹𝑥1 for canard 𝑥1
𝐺 Global frame of reference
I𝑛 𝑛 × 𝑛 identity matrix
𝐽𝑅 Radial mass moment of inertia
𝐽𝐿 Longitudinal mass moment of inertia
J Mass moment of inertia matrix, {𝐽𝑅 𝐽𝐿 𝐽𝐿 } 𝑇 ⋅ I3
K LQR gain matrix
𝐾𝑝 Proportional gain (PID controller)
𝐾𝑖 Integral gain (PID controller)
𝐾𝑑 Derivative gain (PID controller)
𝐿𝑑𝑎𝑚𝑝 Damping moment about roll axis
𝐿𝑐𝑎𝑛 Moment induced by the canards about the roll axis
𝐿𝛼 Canard roll moment derivative w.r.t angle of attack on canard
𝐿⁄ Lift to drag ratio
𝐷
𝑀𝑐𝑜𝑟𝑟 Corrective moment about pitch axis
𝑀𝑑𝑎𝑚𝑝 Damping moment about pitch axis
𝑀𝑐𝑎𝑛 Moment induced by the canards about the pitch axis
𝑀𝛼 Canard pitch moment derivative w.r.t angle of attack on canard

viii
 𝑚̇ Engine exhaust gas mass flow rate
𝐌 Resultant moment vector
𝑁𝑐𝑜𝑟𝑟 Corrective moment about yaw axis
𝑁𝑑𝑎𝑚𝑝 Damping moment about yaw axis
𝑁𝑐𝑎𝑛 Moment induced by the canards about the yaw axis
𝑁𝛼 Canard yaw moment derivative w.r.t angle of attack on canard
𝑁 Number of fins
𝑁 Chernoff bound number of samples
𝜙 Roll angle
𝜃 Pitch angle
𝜓 Yaw angle
𝑝 Roll rate
𝑞 Pitch rate
𝑟 Yaw rate
𝜌 Air density
𝑞̅ Dynamic pressure
Q State-cost weighted matrix
R Input-cost weighted matrix
∑𝑐𝑖 𝜉𝑖2 Δ𝜉𝑖 Roll damping fin planform sum term (Niskanen, 2009)
𝑆 Canard planform area
𝐮 Input vector, {𝑥1 𝑥2 𝑦1 𝑦2 } 𝑇
𝑣 Velocity
𝛚 Attitude rate vector, {𝑝 𝑞 𝑟} 𝑇
𝛚̇ Attitude acceleration vector, {𝑝̇ 𝑞̇ 𝑟̇ } 𝑇
𝑥𝐵 , 𝑦𝐵 , 𝑧𝐵 𝑥, 𝑦, and 𝑧 coordinates in the body frame of reference
𝑥𝐺 , 𝑦𝐺 , 𝑧𝐺 𝑥, 𝑦, and 𝑧 coordinates in the global frame of reference
𝑥1 , 𝑥2 Pitch axis canard deflections
𝑦1 , 𝑦2 Yaw axis canard deflections
𝑋𝐶𝑃 Distance from nose cone tip to centre of pressure
𝑋𝐶𝐺 Distance from nose cone tip to centre of gravity
𝑋̅𝐶𝑃 Distance between CP and CG, 𝑋𝐶𝑃 − 𝑋𝐶𝐺
𝑋𝑇 Distance from nose cone tip to nozzle throat
𝑋̅𝑇 Distance between CG and nozzle throat, 𝑋𝑇 − 𝑋𝐶𝐺
𝑋𝐶𝑃𝑐𝑎𝑛 Distance from nose cone tip to canard CP
𝑋̅𝐶𝑃𝑐𝑎𝑛 Distance between CG and canard CP, 𝑋𝐶𝑃𝑐𝑎𝑛 − 𝑋𝐶𝐺
𝐱 State vector, {𝜙 𝜃 𝜓 𝑝 𝑞 𝑟} 𝑇
𝐱̇ State derivative vector, {𝑝 𝑞 𝑟 𝑝̇ 𝑞̇ 𝑟̇ } 𝑇
A, B, C, D, E, F State space matrices

ix
 𝑌𝐶𝑃𝑐𝑎𝑛 Distance from body frame origin to canard CP in 𝑦𝐵 direction
𝑌̅𝐶𝑃𝑐𝑎𝑛 Distance between CG and canard CP in 𝑦𝐵 direction, 𝑌𝐶𝑃𝑐𝑎𝑛 − 𝑌𝐶𝐺
𝑍𝐶𝑃𝑐𝑎𝑛 Distance from body frame origin to canard CP in 𝑧𝐵 direction
𝑍̅𝐶𝑃𝑐𝑎𝑛 Distance between CG and canard CP in 𝑧𝐵 direction, 𝑍𝐶𝑃𝑐𝑎𝑛 − 𝑍𝐶𝑃

Acronyms

CFD Computational fluid dynamics

CG Centre of gravity
CP Centre of pressure
LQR Linear quadratic regulator (controller)
LURA Leeds University Rocketry Association
MIMO Multi-input, multi-output
PID Proportional-integral-derivative (controller)
RAM Random access memory
SISO Single-input, single-output
SM Static margin

x
 1 Introduction

1.1 Project Aptos

This report forms part of the wider Project Aptos. Conducted as a five-person master’s
team project at the University of Leeds, Project Aptos aims to develop an active control
system to improve the dynamic stability of a sounding rocket and ensure it flies a more
vertical trajectory.

This project was run in conjunction with the Leeds University Rocketry Association
(LURA) and was sponsored by Collins Aerospace. LURA aims to use the findings of this
project to further develop the team’s understanding of active control and will use future
iterations of the Aptos system to ensure their rockets fly vertically, and therefore reach
higher altitudes, in the pursuit of the UK amateur rocketry altitude record (UKRA, 2021).

Following consideration of different control systems, it was decided to use a set of four
canards near the nose of the rocket controlled by servomotors to steer the rocket during
ascent. Details of the full system, including the Pathfinder rocket on which the control
system was implemented, are discussed by Youds et al. (2023).

1.2 Individual Aims and Objectives

This report discusses the design, analysis, and testing of the Aptos control algorithm.
The aim of this project was to develop a closed loop controller capable of making the
rocket fly a more vertical trajectory by controlling the deflections of the four canards.

This aim consisted of several objectives:

1 Create a mathematical model of the system dynamics from first principles.

2 Implement this model in Simulink in order to predict the behaviour of the
system.
3 Design a controller that can deflect the canards in order to maintain more
vertical flight.

1.3 Report Structure

These objectives are addressed in order through this report. First, existing literature on
rocket active control systems is reviewed in Chapter 2. A mathematical model of the
aerodynamic behaviour of the rocket is derived in Chapter 3, and this is then used to
design a closed loop controller in Chapter 4. The position control of the servomotors is
also discussed here. Chapter 4 concludes with analysis of the system’s performance
and robustness to expected variations in the rocket’s geometry.

1
2 Literature Review and Theoretical Background

2.1 Sounding Rocket Stability

A rocket is considered passively stable if its centre of pressure (CP) is further aft than its
centre of gravity (CG), such that the aerodynamic forces acting on the airframe keep the
rocket aligned with the incoming airflow (NAROM Andøya Space Centre, 2018). It is
generally accepted that a rocket must maintain a static margin, 𝑆𝑀 > 1, to be passively
stable, where the static margin is given by the distance between the CP and the CG
divided by the rocket’s body diameter.

Passive stability is not always desirable, especially when a rocket aims to reach high
altitudes. As identified by Uselton (1971) and Ujjin et al. (2021), passively stable rockets
‘weathercock’ in the presence of a crosswind. Whilst not flying vertically, the horizontal
component of the thrust vector is effectively wasted, limiting altitude and increasing the
possible landing radius (Lorenz and Bierig, 2013). What is more desirable for such
rockets is an active system, which maintains dynamic stability and allows passively
unstable rockets to fly pre-determined trajectories.

Active aerodynamic control surfaces can be used to steer rockets in a similar fashion to
conventional aircraft. Their use is widely documented on missile and artillery systems,
for example by Guo et al. (2016), Uselton (1971) and Seaberg and Smith (1951), where
they are generally used for precise impacts. Despite the different use case, some control
laws designed for weapon systems can be applied to sounding rockets.

Mracek and Ridgely (2006), for example, discuss at depth the autopilot system for a
canard-controlled missile. Their system has great similarities with the system developed
in this project, but the missile uses two sets of control surfaces: canards and rear fins.

Similar research exists for non-military rockets, although it generally focuses on larger-
scale vehicles than Pathfinder. Choi and Bang (2000) developed an adaptive controller
specifically for a slender sounding rocket, where the body-bending vibration of the
airframe is considered and compensated for by the controller. Despite the similarities
between Pathfinder and the KSR-II rocket examined in the paper, this level of structural
analysis is not required for this project, where the fibreglass airframe provides good
rigidity relative to the low thrust and aerodynamic loads.

Braswell et al. (2017) discuss the development of a canard control system as a payload
for the NASA Student Launch Competition. A series of flight tests were performed to
tune the controller, though a loose fin on the competition flight prevented the system

2
from performing as expected. The scale of the rocket itself is not dissimilar to Pathfinder
and the budget and timescales align well with this project, though their system focusses
solely on controlling roll, in contrast to the three-axis control needed for Aptos.

2.2 Closed Loop Control Systems

These papers provide a good insight into different control law options. Since Braswell et
al. (2017) aim to control roll attitude only, the system is single-input, single-output
(SISO), and thus a PID controller is used. This controller is simple to setup and tuning it
does not require a linear model of the system, though the effects of the gain values
themselves are not always intuitive for more complex nonlinear systems.

Another option is an LQR, as implemented by Mracek and Ridgely (2006). LQR

controllers are inherently robust, handle multi-input, mutli-output (MIMO) systems very
well, and are intuitive to tune. Actuating both the canards and rear fins introduces
additional complexity in Mracek and Ridgely’s work, though their strategy for tuning and
analysing the system’s performance was useful for this project. Aside from rockets,
discussion of various control systems for quadrotor drones (Vigneswaran and Kp, 2019)
and satellites (Narkiewicz et al., 2020) also provided helpful research for this project.

2.3 Robustness Analysis

Variations in the rocket’s geometry and mass properties occur during the manufacturing
process and therefore in reality, the system dynamics will differ from the model that the
controller has been designed with. It is therefore important to verify that a controller is
robust to these variations in terms of stability and performance.

Robustness and stability can be investigated by examining the worst perturbation from
the nominal state space in the form of 𝜇 analysis (Morton and McAfoos, 1985), or by
calculating the disk margins of the system (Seiler et al., 2020). These methods provide
a deterministic solution to the robustness problem, but do not account for the additional
complexity of nonlinear behaviour, which is inherent to any aerodynamic system.

Alternatively, a more stochastic methodology can be used in the form of the Monte Carlo
approach (Ray and Stengel, 1993). Generally, this method is computationally expensive
and risks missing a parameter variation combination that can make the system unstable.
However, work by Postlethwaite et al. (2009) and Tempo and Dabenne (2004) applies
the Chernoff Bound to calculate the lower bound for the number of samples which must
be used to gain a given accuracy and confidence in a system’s robustness. This proved
extremely useful for the analysis conducted in Chapter 4.

3
 3 Mathematical Modelling

Before a controller can be designed, a mathematical model of the system must be

derived. This model must accurately describe the behaviour of the system dynamics
within the expected operating conditions without being overly computationally
expensive.

3.1 Frames of Reference

Two frames of reference are used here. The body frame, 𝐵, is fixed to the rocket during
flight. It is a right-handed coordinate system with origin at the rocket’s CG, and the 𝑥𝐵
axis parallel to the rocket’s longitudinal axis, pointing upwards towards the nose. The 𝑦𝐵
and 𝑧𝐵 axes are then orthogonal and lie on the rocket’s cross section.

The global frame, 𝐺, is a right-handed coordinate system with its origin fixed to the
ground at the location of the launch pad. The 𝑥𝐺 axis points directly upwards,
perpendicular to the ground and the 𝑦𝐺 and 𝑧𝐺 axes are orthogonal to it parallel to the
ground. The flight is low-altitude and lasts less than a minute, and therefore the Earth’s
rotation has been neglected and the ground has been modelled as a flat plane.

At the time of ignition, the body and global frames are aligned, but do not remain so
during flight. The rocket’s roll, pitch and yaw angles, 𝜙, 𝜃, and 𝜓, are then defined as the
Euler angles that compose the rotation, R, such that any vector, 𝐚 ∈ ℝ3 , is transformed
from the global frame to the body frame, as in (3.1).

𝐚𝐵 = R(𝜙, 𝜃, 𝜓) 𝐚𝐺 (3.1)

3.2 Pitch and Yaw Moments

It has been assumed that there are three moments acting on the rocket’s pitch and yaw
axes during flight.

𝑀𝑐𝑜𝑟𝑟 Corrective moment induced by the lift acting on the rocket’s passive fins.
𝑀𝑑𝑎𝑚𝑝 Damping moment induced by the drag on the airframe as it rotates.
𝑀𝑐𝑎𝑛 Canard moment induced by the lift of the deflected canards.

3.2.1 Pitch / Yaw Corrective Moment

Providing that the rocket’s CP is further aft than its CG, there will be a corrective moment
acting to align the rocket with the incoming air flow. In the presence of a crosswind, this
moment will reduce the rocket’s angle of attack. Mandell et al. (1973) state that the pitch
and yaw corrective moments, 𝑀𝑐𝑜𝑟𝑟 and 𝑁𝑐𝑜𝑟𝑟 , acting on the rocket are proportional to
the angle of attack as in (3.2).

4
 𝑀𝑐𝑜𝑟𝑟 = 𝐶1 𝛼 , 𝑁𝑐𝑜𝑟𝑟 = 𝐶1 𝛽 (3.2)

Where 𝐶1 is the corrective moment coefficient, given by (3.3) and derived by Mandell et
al. (1973).

𝐶1 = 𝑞̅𝐴𝑟𝑒𝑓 𝐶𝑁𝛼 𝑋̅𝐶𝑃 (3.3)

1
Where 𝑞̅ is the dynamic pressure of the airflow, 𝑞̅ = 2 𝜌𝑣 2 , 𝐴𝑟𝑒𝑓 is the body tube’s cross-

sectional area, and 𝑋̅𝐶𝑃 is the moment arm between the CP and the CG:
𝑋̅𝐶𝑃 = 𝑋𝐶𝑃 − 𝑋𝐶𝐺 . 𝐶𝑁𝛼 is the rocket’s normal force coefficient derivative with respect to
angle of attack. The methods used to calculate this parameter were derived by Niskanen
(2009) and Youds (2022) and are included in Appendix A.

3.2.2 Pitch / Yaw Damping Moment

Mandell et al. (1973) also derive the rocket’s pitch and yaw damping moments, and as
expected, these are proportional to the pitch rate, 𝑞, and yaw rate, 𝑟, respectively.

𝑀𝑑𝑎𝑚𝑝 = 𝐶2 𝑞 , 𝑁𝑑𝑎𝑚𝑝 = 𝐶2 𝑟 (3.4)

Where 𝐶2 is the damping moment coefficient, derived by Mandell et al. (1973).

𝑞̅ 2
𝐶2 = 𝐴 [∑𝐶𝑁𝛼 𝑋̅𝐶𝑃 𝑖 ] + 𝑚̇ 𝑋̅𝑇2 (3.5)
𝑣 𝑟𝑒𝑓 𝑖

Where 𝐶𝑁𝛼 is the normal force coefficient derivative of the 𝑖 th fuselage component, (e.g.,
𝑖

nose cone, body tube, or fins), and 𝑋̅𝐶𝑃𝑖 is the moment arm between the rocket’s CG
and the CP of the 𝑖 th component: 𝑋̅𝐶𝑃𝑖 = 𝑋𝐶𝑃𝑖 − 𝑋𝐶𝐺 . The first term here represents the
aerodynamic damping due to drag as the rocket rotates, and the second represents the
jet damping.

Jet damping is an additional damping torque induced by the jet of exhaust gases leaving
the engine. This fast-moving gas resists changes in rotational momentum and helps to
keep the rocket flying straight. 𝑚̇ is the mass flow rate of the exhaust gas, and 𝑋̅𝑇 is the
moment arm between the rocket’s CG and the engine throat.

3.2.3 Pitch / Yaw Canard Moment

When the canards are deflected, the lift force they generate induces pitch and yaw
moments on the rocket. Figure 5 in Appendix D shows the canard configuration.

Two canards on opposite sides of the rocket control pitch motion about the 𝑦𝐵 axis. The
canard moment, 𝑀𝑐𝑎𝑛 , is equal to the distance between the rocket’s CG and the CP of
the canards multiplied by the horizontal component of the canard lift force. Drag can be

5
assumed to be negligible here since the canard 𝐿⁄𝐷 ratio is approximately 120 (Daney
de Marcillac, 2023). Therefore, the canard pitch moment is given by (3.6).

𝑀𝑐𝑎𝑛 = 𝑞̅𝑆𝑋̅𝐶𝑃𝑐𝑎𝑛 (𝐶𝐿1 cos 𝑥1 + 𝐶𝐿2 cos 𝑥2 ) (3.6)

The lift coefficient itself is proportional to angle of attack. The angle of attack on the
canard depends on the rocket’s angle of attack, 𝛼, and the canard deflection, 𝑥𝑖 . The
sign convention for the canard deflections is important here. Since the body frame origin
lies at the CG, one canard lies on the 𝑦𝐵+ axis – canard 1 – and the other on the 𝑦𝐵− axis
– canard 2. The deflections of these canards are 𝑥1 and 𝑥2 respectively.

A positive 𝑥1 deflection generates a positive pitching moment, as per the righthand screw
rule. The opposite is true for canard 2; a positive 𝑥2 deflection results in a negative
pitching moment. Therefore, if 𝑥1 = 𝑥2 , there will be no pitching motion, only rolling
motion. Consequently, the AoA on canard 1 is the difference between 𝛼 and 𝑥1 , whereas
for canard 2, it is the sum of 𝛼 and 𝑥2 . Applying this to (3.6) gives (3.7).

𝐶𝐿𝑖 = 𝐶𝐿𝛼 𝛼𝑐𝑎𝑛 𝑖

⇒ 𝑀𝑐𝑎𝑛 = 𝑞̅𝑆𝐶𝐿𝛼 𝑋̅𝐶𝑃𝑐𝑎𝑛 [(𝛼 − 𝑥1 ) cos 𝑥1 + (𝛼 + 𝑥2 ) cos 𝑥2 ] (3.7)

For simplicity of notation, the pitch moment derivative with respect to canard angle of
attack can be defined as 𝑀𝛼 = 𝑞̅𝑆𝐶𝐿𝛼 𝑋̅𝐶𝑃𝑐𝑎𝑛 . Applying this substitution to (3.7) gives (3.8).

𝑀𝑐𝑎𝑛 = 𝑀𝛼 [(𝛼 − 𝑥1 ) cos 𝑥1 + (𝛼 + 𝑥2 ) cos 𝑥2 ] (3.8)

This is also true for the yaw canard moment, 𝑁𝑐𝑎𝑛 .

𝑁𝑐𝑎𝑛 = 𝑁𝛼 [(𝛽 − 𝑦1 ) cos 𝑦1 + (𝛽 + 𝑦2 ) cos 𝑦2 ] (3.9)

3.3 Roll Moments

Two moments act on the rocket’s roll axis. Since the Pathfinder fins are mounted on
bearings and can spin freely, it has been assumed that the only roll forcing moment is
due to the canards.

𝐿𝑑𝑎𝑚𝑝 Damping moment induced by the drag on the airframe as it rotates.

𝐿𝑐𝑎𝑛 Canard moment induced by the lift of the deflected canards.

3.3.1 Roll Damping Moment

Youds (2022) expands Niskanen’s (2009) work to show that the roll damping moment is
proportional to the rocket’s roll rate, 𝑝, as in (3.16).

𝐿𝑑𝑎𝑚𝑝 = 𝐶𝑑𝑝 𝑝 (3.10)

Where 𝐶𝑑𝑝 is the roll damping coefficient with respect to roll rate (Youds, 2022).

6
 𝑞̅
𝐶𝑑𝑝 = 𝐴 𝑑 𝑁𝐶 (∑𝑐𝑖 𝜉𝑖2 Δ𝜉𝑖 ) (3.11)
𝑣 𝑟𝑒𝑓 𝑏𝑜𝑑𝑦 𝑁𝛼0
Where 𝑑𝑏𝑜𝑑𝑦 is the body tube diameter and 𝐶𝑁𝛼0 is the normal force coefficient of the

fin’s 2D aerofoil. The final sum term is constant for a specific fin planform shape and is
given by Youds (2022) and Niskanen (2009).

3.3.2 Roll Canard Moment

The canard roll moment is defined by (3.12). Figure 6 in Appendix D shows the directions
of these forces. Again, the drag on the canards has been neglected.

𝐿𝑐𝑎𝑛 = 𝑌̅𝐶𝑃𝑐𝑎𝑛 𝐹𝑥1 − 𝑌̅𝐶𝑃𝑐𝑎𝑛 𝐹𝑥2 + 𝑍̅𝐶𝑃𝑐𝑎𝑛 𝐹𝑦1 − 𝑍̅𝐶𝑃𝑐𝑎𝑛 𝐹𝑦2 (3.12)
Where 𝑌̅𝐶𝑃𝑐𝑎𝑛 and 𝑍̅𝐶𝑃𝑐𝑎𝑛 are the moment arms between the CG and the CP of the
canards in the 𝑦𝐵 and 𝑧𝐵 axes respectively. It is assumed that the CG lies on the body
frame origin, and therefore 𝑌̅𝐶𝑃𝑐𝑎𝑛 = 𝑌𝐶𝑃𝑐𝑎𝑛 and 𝑍̅𝐶𝑃𝑐𝑎𝑛 = 𝑍𝐶𝑃𝑐𝑎𝑛 . Moreover, the canards
are placed symmetrically about the 𝑥 axis, and thus 𝑌𝐶𝑃𝑐𝑎𝑛 = 𝑍𝐶𝑃𝑐𝑎𝑛 .

𝐹𝑥1 , 𝐹𝑥2 , 𝐹𝑦1 , and 𝐹𝑦2 are the components of the lift forces acting in the 𝑦𝐵 − 𝑧𝐵 plane.
Since 𝐹𝑖 = 𝑞̅𝑆𝐶𝐿𝛼 𝛼𝑐𝑎𝑛 𝑖 cos 𝑥𝑖 , and using the angles of attack derived in Section 3.2.3, the
total moment can be expressed in the form shown in (3.13).

𝐿𝑐𝑎𝑛 = 𝐿𝛼 [(𝛼 − 𝑥1 ) cos 𝑥1 − (𝛼 + 𝑥2 ) cos 𝑥2 + (𝛽 − 𝑦1 ) cos 𝑦1 − (𝛽 − 𝑦2 ) cos 𝑦2 ] (3.13)

Where 𝐿𝛼 = 𝑞̅𝑆𝐶𝐿𝛼 𝑌𝐶𝑃𝑐𝑎𝑛 = 𝑞̅𝑆𝐶𝐿𝛼 𝑍𝐶𝑃𝑐𝑎𝑛 is the derivative of the canard roll moment with
respect to the angle of attack of the respective canard.

3.4 Equation of Motion Formulation

Euler’s Rotational Equation of Motion relates the rocket’s rotational accelerations, 𝛚̇,
and rates, 𝛚, to the moments acting on it during flight, as shown in (3.14).

J𝛚̇ + 𝛚 × (J𝛚) = 𝐌 (3.14)

Where J is the rocket’s moment of inertia matrix. It can be assumed that the rocket is
axisymmetric about the 𝑥𝐵 axis, and therefore J is given by (3.15).

J = {𝐽𝑅 𝐽𝐿 𝐽𝐿 } 𝑇 ⋅ I3 (3.15)

Where 𝐽𝑅 and 𝐽𝐿 are the rocket’s radial and longitudinal mass moments of inertia
respectively and I𝑛 is the 𝑛 × 𝑛 identity matrix.

The resultant moment vector, 𝐌, is given by (3.16).

𝐿𝑐𝑎𝑛 − 𝐿𝑑𝑎𝑚𝑝
𝐌 = [ 𝑐𝑎𝑛 − 𝑀𝑐𝑜𝑟𝑟 − 𝑀𝑑𝑎𝑚𝑝 ]
𝑀
𝑁𝑐𝑎𝑛 − 𝑁𝑐𝑜𝑟𝑟 − 𝑁𝑑𝑎𝑚𝑝

7
 𝐿𝛼 [(𝛼 − 𝑥1 ) cos 𝑥1 − (𝛼 + 𝑥2 ) cos 𝑥2 + (𝛽 − 𝑦1 ) cos 𝑦1 − (𝛽 + 𝑦2 ) cos 𝑦2] − 𝐶𝑑𝑝 𝑝
⇒𝐌=[ 𝑀𝛼 [(𝛼 − 𝑥1 ) cos 𝑥1 + (𝛼 + 𝑥2 ) cos 𝑥2 ] − 𝐶1 𝛼 − 𝐶2 𝑞 ] (3.16)
𝑁𝛼 [(𝛽 ) (𝛽 ) ]
− 𝑦1 cos 𝑦1 + + 𝑦2 cos 𝑦2 − 𝐶1 𝛽 − 𝐶2 𝑟

Equations (3.15) and (3.16) were substituted into (3.14) to give the rotational equations
of motion for the rocket. The final equation has been omitted here for the sake of brevity.

3.5 Nonlinear Modelling

To run simulations of the rocket’s motion and tune the controller, the nonlinear dynamics
equation (3.14) was implemented in Simulink. Aerodynamic properties, such as the
normal force coefficient derivatives and centres of pressure were calculated directly from
Pathfinder’s geometry using the equations discussed in Appendix A and the work of
Youds (2022) and Niskanen (2009).

Mass properties, such as the mass moments of inertia, centre of gravity location, and
exhaust gas mass flow rate were exported from OpenRocket. Lookup tables were
implemented in Simulink to model the variation of these parameters throughout ascent.

The vertical velocity over time was also exported from OpenRocket and a range of
profiles were tested to ensure the controller was not overly sensitive to velocity variations
as this is something that was expected to be difficult to perfectly predict for the real flight.

The complexity of the modelling was increased over time. First, the system response
was analysed at a single operating point, for example a constant velocity and altitude.
This was used to ensure the system had adequate controllability before progressing.

A time variant model was then created, where the operating condition was itself a
function of time, to better simulate a real flight. To start with, this model assumed a
crosswind of constant speed and direction, acting only in the 𝑦𝐺 − 𝑧𝐺 plane.

Later, this was improved, using Dryden’s wind turbulence model in Simulink
(MathWorks, n.d.). This model varied both the windspeed and direction based on the
rocket’s altitude and attitude. This was used to ensure the rocket could perform three
axis control effectively and respond adequately to sudden gusts of wind. Full analysis of
these results has not been included in the main body of the report due to space
limitations, but graphs of the system performance for these models is in Appendix E.

8
 4 Closed Loop Controller Design

4.1 States, Inputs, and Disturbances

To fly vertically, the controller works to minimise the rocket’s roll, pitch, and yaw angles
and rates. These six variables form the state vector, 𝐱 ∈ ℝ𝟔 , shown in (4.1), where 𝜙, 𝜃,
and 𝜓 are the roll, pitch and yaw angles respectively, and 𝑝, 𝑞, and 𝑟 are their rates.

𝐱 = {𝜙 𝜃 𝜓 𝑝 𝑞 𝑟} 𝑇 (4.1)

These four angular canard deflections make up the input vector, 𝐮 ∈ ℝ𝟒 , shown in (4.2).
The notation for the canard naming described in Chapter 3 is used again here.

𝐮 = {𝑥1 𝑥2 𝑦1 𝑦2 } 𝑇 (4.2)

As shown in Chapter 3, the angle of attack and sideslip angle, 𝛼 and 𝛽, affect the lift
generated by the canards, as well as the pitch and yaw corrective moments. Measuring
these angles would require two angle of attack vanes as seen on aircraft, placed
perpendicular to each other on the rocket’s airframe (Gracey, 1958). Those used on
aircraft are too large and heavy for use on Pathfinder, and therefore custom sensors
would need to be designed, manufactured and tested, adding significant complexity.

Instead, 𝛼 and 𝛽 are modelled as disturbances and are not measured. The disturbance
vector, 𝐝 ∈ ℝ𝟐 , is therefore of the form shown in (4.3).

𝐝 = {𝛼 𝛽} 𝑇 (4.3)

4.2 State Space Formulation

The rocket’s rotational equations derived in Chapter 3 are inherently nonlinear. A linear
representation can be formed, known as the state space formulation, where the system
dynamics are described by equation (4.4).
𝐱̇ = A𝐱 + B𝐮 + E𝐝
(4.4)
𝐲 = C𝐱 + D𝐮 + F𝐝

Where A, B, C, D, E, and F are matrices that relate the system states, inputs, and
disturbances to the state derivative, 𝐱̇ , and output, 𝐲.

To obtain this formulation, the dynamics equations derived in Section 3.4 were
rearranged in terms of the components of 𝐱, 𝐮, and 𝐝. These equations were then
linearised around the steady-state operating point 𝐱̅ = 0, using the first-order Taylor
series and small angle approximations. Full derivation of the state space equation is
presented as Appendix B. The final state space representation is shown as (4.5).

9
 0 0 0 1 0 0 0 0 0 0
0 0 0 0 1 0 0 0
0 0 0 0
0 0 0 0 0 1 0 0
0 0 0 0
𝐶𝑑𝑝 𝐿𝛼 𝐿𝛼 𝐿𝛼 𝐿𝛼 0 0
0 0 0 − 0 0 − − − − 0 0
𝐱̇ = 𝐽𝑅 𝐱+ 𝐽𝑅 𝐽𝑅 𝐽𝑅 𝐽𝑅 𝐮 + 𝑀𝛼 𝐶1 𝐝(4.5)
𝐶2 𝑀𝛼 𝑀𝛼 2 − 0
0 0 0 0 − 0 − + 0 0 𝐽𝐿 𝐽𝐿
𝐽𝐿 𝐽𝐿 𝐽𝐿 𝑁𝛼 𝐶1
𝐶2 𝑁𝛼 𝑁𝛼 0 2 −
0 0 0 0 0 − 0 0 − + [ 𝐽𝐿 𝐽𝐿 ]
[ 𝐽𝐿 ] [ 𝐽𝐿 𝐽𝐿 ]
𝐲 = I6 𝐱

4.3 LQR Design

Once the state space representation was formed, the LQR controller was tuned. Figure
1 shows the layout of the closed loop control system.

Figure 1 - Closed loop control diagram of the Aptos control system.

The LQR controller acts as a 4×6 pure gain matrix, K, applied to the state error in order
to derive the desired control inputs, 𝐮. These gains are calculated by minimising the cost
function (4.6), where 𝐮 = −K 𝐱 and Q and R are the state-cost and input-cost weighted
matrices (Yang, 2012; Okyere et al., 2019).
∞
𝐽 = ∫ (𝐱 T Q𝐱 + 𝐮T R𝐮) 𝑑𝑡 (4.6)
0

From inspection, this is an energy minimisation problem, where 𝐱 T Q𝐱 represents the

energy of the rocket due to its current state, and 𝐮T R𝐮 represents the energy associated
with the actuation of the inputs.

In practice, to tune the controller, the MATLAB lqr(A, B, Q, R) function was used,
where A and B represent the state space matrices from (4.5). The Q and R matrices are
varied to achieve the desired performance. Since each canard is actuated by an identical
servomotor, the cost was considered equal across all inputs and the input-cost weighted
matrix was of the form shown in (4.6).

10
 R = 𝑟 ⋅ I4 (4.6)

Where 𝑟 is the actuator penalisation weight.

This was not the case for the state-cost weighted matrix, where individual weights were
used to prioritise the reduction of some state errors more than others, as in (4.7).

Q = {𝑞𝜙 𝑞𝜃 𝑞𝜓 𝑞𝑝 𝑞𝑞 𝑞𝑟 } 𝑇 ⋅ I6 (4.7)

Where 𝑞𝑖 are the weightings for each state.

For the first few flight tests, a priority was to ensure that the canard motion was as
smooth as possible, avoiding sudden changes in deflection. This helps prevent the
canards from undergoing dynamic stall and other nonlinear aerodynamic effects.

To achieve this, minimising the pitch and yaw angles was prioritised over minimising the
pitch and yaw rates. Thus, 𝑞𝜃 > 𝑞𝑞 and 𝑞𝜓 > 𝑞𝑟 . As a result, the response was slightly
slower acting and less damped, however the canard deflections were smoother.

During flight, minimising roll rate is significantly more important than roll angle. If the
rocket begins to spin at a fast rate, a small pitch or yaw perturbation will lead to a
significant gyroscopic precession known as coning (Youds, 2022). Pitch and yaw
controllability will be significantly reduced to the point that instability could occur.

When the rocket rolls, it’s important that the controller reduces the roll rate, but does not
attempt to return the rocket to its original roll angle. To this end, the ideal roll angle
weighting, 𝑞𝜙 , is zero. However, one condition of the LQR controller is that both the Q
and R matrices must be positive definite (Yang, 2012) and therefore a very small value
is used such that 𝑞𝑝 ≫ 𝑞𝜙 .

The specific values of the 𝑞 and 𝑟 weighting factors were determined through trial-and-
error until the Simulink simulations indicated that the system met the performance
requirements. Bryson’s Rule (Okyere et al., 2019) was investigated as a method to
achieve this algorithmically, although the gains this produced were very high and it was
later found that comparable performance could be achieved with significantly lower
gains. The final weightings are shown in (4.8) and (4.9).

Q = {1 × 10−9 4 4 4 1 1 } 𝑇 ⋅ I6 (4.8)

R = 5 ⋅ I4 (4.9)

4.4 Gain Scheduling

Every aerodynamic property of the system varies with velocity, and most are proportional
to the velocity squared. Intuitively, the controller will require significantly less actuator
effort when the rocket is flying at higher speeds than it does when at low speeds.
11
Gain scheduling was implemented to ensure that the controller could account for these
velocity changes as the rocket accelerates. Optimum gain-sets were calculated at
different velocities from 30 m/s up to 150 m/s (about 1.5x the expected maximum speed).
The onboard flight computer then hands off between gain-sets as the differentiated
barometer data indicates that the vertical velocity has changed. Gain scheduling
significantly improved performance across the full velocity profile with very little memory
cost on the flight computers.

4.5 Servomotor Position Control

The analysis up to this point has assumed that the LQR controller can command a
certain canard deflection, and this deflection is achieved instantly. In reality, both time
delay and position error are present in the servomotors, and this affects the system’s
performance.

The Herkulex DRS-0101 motors used in the Aptos module were smart servos and
featured an in-built potentiometer and PID controller for position control. The PID gains
could be written directly to the servos’ RAM, meaning that different combinations could
be tested quickly and easily.

From requirements APT-REQ-1-5-1 and APT-REQ-1-5-6, the servomotors must

respond to a 10° step input with a settling time less than 0.1 s and an accuracy of ±0.5°
with no overshoot outside of these bounds. In addition, a pointing accuracy of 1±0.1° is
needed to ensure the servos have the fine control required for flight (APT-REQ-1-5-5).

This was achieved through a trial-and-error testing process. One of the servos was
connected to a computer via an Arduino Uno development board and a step input was
commanded. The angular position of the servo was read from its potentiometer, and the
PID gains were iteratively tuned until performance met the requirements.

Figure 2 shows the response to both a 1° and 10° step input. The PID gains used are
displayed in Table 1.

Figure 2 - Servomotor position step response graphs with the default and updated PID gain sets for steps
of 1° and 10°.
12
 Table 1 - PID gains used for the servomotor position controller.

𝑲𝒑 𝑲𝒊 𝑲𝒅
Default 440 0 8000
Updated 550 1000 9000

The default gains provided a good response for the 10° step, however the performance
for the 1° step was poor. The servomotors operate in increments of 0.325°, and the gains
were not high enough to overcome the internal friction needed to step through the final
increment, thus the shaft gets stuck at 0.65°. Similar plateaus can be seen in the updated
gain response, however the higher 𝐾𝑝 and 𝐾𝑖 gains forced the controller to overcome
the frictional force.

The updated response for the 10° step was very good. There was a small amount of
overshoot, but this was within acceptable bounds. Increasing the 𝐾𝑑 gain further reduced
this overshoot but slowed the response for the 1° step too much.

Using the MATLAB System Identification App and the step data shown in Figure 2, a
transfer function was fit, which represented the servomotor position control dynamics.
This was used in Simulink to model the delay and position error associated with actuating
the canards.

4.6 Robustness Analysis

Requirement APT-REQ-1-1-1 states that the controller must be robust to expected

deviations from the geometry of the Pathfinder rocket with 95% accuracy and
confidence. Specifically, the nominal values and expected variability of each of these
parameters are provided in Appendix C.

The Chernoff Bound can be used to calculate the number of samples needed to gain the
required accuracy and confidence, as in (4.8) (Postlethwaite et al., 2009).

log(2⁄𝛿 )
𝑁≥ (4.8)
2𝜖 2
Applying 𝜖 = 𝛿 = 0.05 gives a minimum sample number of 737.8. In total, 740
simulations were run. Figure 3 shows how the system is robust both in terms of
performance and stability across the entire expected parameter variation. This model
assumes a constant windspeed of 5 m/s acting only in the 𝑧𝐺+ direction.

This is expected as the LQR is an inherently robust controller, however it is impressive

that the system can handle up to ±50% variation in the canard 𝐶𝐿𝛼 value. This fact
prevented the need for any high fidelity aerodynamic modelling of the canards.

13
 Figure 3 - Results from the control system robustness analysis.

4.7 Controller Performance

Figure 4 displays how the system handles a 5 m/s constant crosswind. In this simulation,
the wind vector acted solely in the 𝑧𝐺+ direction and therefore induced pitch motion only.
For comparison, the free response without active control enabled is also plotted on
Figure 4.

Figure 4 - System performance with 5 m/s constant crosswind.

The control algorithm was set to only activate whilst the vertical velocity, 𝑣 ≥ 30 m/s.
This prevented it from activating whilst the rocket was still on the launch rail and
minimises actuator energy cost towards apogee when the velocity is too low for the
rocket to be controllable.

The Aptos system reduced the peak pitch angle by about 2.7° and significantly damped
the oscillations following rail departure. The steady state pitch angle during ascent is
very similar for both cases, as is the increase in pitch as the rocket decelerates towards
apogee. This system provides a good first proof-of-concept, and performance can be
significantly improved by reducing the input-cost weighting matrix, R.
14
 5 Conclusion

5.1 Achievements and Impact

This project aimed to develop the active control algorithm for a canard-controlled
sounding rocket stability augmentation system. First, the work of Mandell et al. (1973),
Niskanen (2009), and Youds (2022) was expanded to derive a three degree of freedom
dynamic model of a canard-controlled sounding rocket. This model was linearised using
first order Taylor series techniques and small angle approximations, and an LQR
controller was tuned to minimise the rocket’s attitude and attitude rates throughout
ascent. A separate PID controller was tuned to improve the response of the servomotors
that control the canards, and system identification techniques were used to model the
actuator’s behaviour in Simulink. Finally, the system’s performance was verified and its
robustness was analysed using stochastic methods derived by Postlethwaite et al.
(2009) and Tempo and Dabenne (2004).

The performance of the LQR controller in simulations indicated that it met the
requirements set out in the project’s contract performance plan. Moreover, due to the
controller’s inherent robustness, it has been proven to handle significant variations in
rocket and canard geometry with almost no cost to performance or stability. Specifically,
this was verified with both an accuracy and confidence of 95%.

Significant progress has been made towards developing a fully functional, modular
canard control system. Although full flight testing has not yet been possible, the two data
acquisition flights conducted throughout this project were highly successful. The
performance of the control algorithm during both flights was nominal and indicates a high
likelihood of success in the future. More discussion on the control algorithm’s
performance during these flights is included in Youds et al. (2023).

5.2 Future Work

More flight tests are already planned for the remainder of 2023, specifically aiming to fly
several times with the active control system enabled. Future work will make use of the
existing datasets to improve confidence in the system and increase performance further.
For example, adaptive control techniques could be employed, using real flight data to
tune the next iteration of the Aptos control algorithm.

It is also recommended that hardware-in-the-loop testing is conducted using flight data

as further improvements in performance can likely be found in the optimisation of the
embedded software running the algorithm.

15
 6 References
Braswell, V., Arcenaux, H., Strutzenberg, H., Taylor, R. and Lineberry, D.M. 2017. UAH
Student Launch Rocket Project: Roll Control of a Sounding Rocket Through
Aerodynamic Surfaces In: 53rd AIAA/SAE/ASEE Joint Propulsion Conference
[Online]. Atlanta, GA: American Institute of Aeronautics and Astronautics.
[Accessed 2 December 2021]. Available from:
[Link]

Choi, H.-D. and Bang, H. 2000. An adaptive control approach to the attitude control of a
flexible rocket. Control Engineering Practice. 8(9), pp.1003–1010.

Daney de Marcillac, J. 2023. Design and Aerodynamic Analysis of a Canard-Controlled

Sounding Rocket. Leeds, UK: University of Leeds.

Gracey, W. 1958. Summary of Methods of Measuring Angle of Attack on Aircraft.

Available from:
[Link]

Guo, Q., Song, W., Gao, M. and Fang, D. 2016. Advanced Guidance Law Design for
Trajectory-Corrected Rockets with Canards under Single Channel Control.

James Barrowman 1967. The Practical Calculation of the Aerodynamic Characteristics

of Slender Finned Vehicles.[Online] Washington D.C.: The Catholic University of
America. [Accessed 13 December 2021]. Available from:
[Link]

Lorenz, S. and Bierig, A. 2013. Robustness Analysis Related to the Control Design of
the SHEFEX-II Hypersonic Canard Control Experiment In: AIAA Guidance,
Navigation, and Control (GNC) Conference [Online]. Boston, MA: American
Institute of Aeronautics and Astronautics. [Accessed 1 May 2023]. Available
from: [Link]

Mandell, G.K., Caporaso, G.J. and Bengen, W.P. 1973. Topics in advanced model
rocketry. Cambridge, Mass: MIT Press.

MathWorks n.d. Dryden Wind Turbulence Model (Continuous). MATLAB and Simulink
Documentation. [Online]. Available from:
[Link]
.html.

Morton, B.G. and McAfoos, R.M. 1985. A Mu-Test for Robustness Analysis of a Real-
Parameter Variation Problem In: 1985 American Control Conference [Online].
Boston, MA, pp.135–138. Available from:
[Link]

Mracek, C. and Ridgely, D. 2006. Optimal Control Solution for Dual (Tail and Canard)
Controlled Missiles In: AIAA Guidance, Navigation, and Control Conference and
Exhibit [Online]. Keystone, Colorado: American Institute of Aeronautics and
Astronautics. [Accessed 1 May 2023]. Available from:
[Link]

16
Narkiewicz, J., Sochacki, M. and Zakrzewski, B. 2020. Generic Model of a Satellite
Attitude Control System. International Journal of Aerospace Engineering. 2020,
pp.1–17.

NAROM Andøya Space Centre 2018. Aerodynamics and forces acting on the rocket.
NAROM Student Rocket Pre-Study. [Online]. [Accessed 21 June 2021].
Available from: [Link]
theory/rocket-dynamics-2/4-1aerodynamics-and-forces-acting-on-the-rocket/.

Niskanen, S. 2009. Development of an Open Source model rocket simulation

software.[Online] Helsinki: Helsinki University of Technology. [Accessed 28 July
2021]. Available from: [Link]

Okyere, E., Bousbaine, A., Poyi, G.T., Joseph, A.K. and Andrade, J.M. 2019. LQR
controller design for quad‐rotor helicopters. The Journal of Engineering.
2019(17), pp.4003–4007.

Postlethwaite, I., Kim, J. and Bates, D.G. 2009. A geometrical formulation of the μ-lower
bound problem. IET Control Theory & Applications. 3(4), pp.465–472.

Ray, L.R. and Stengel, R.F. 1993. A monte carlo approach to the analysis of control
system robustness. Automatica. 29(1), pp.229–236.

Seaberg, E. and Smith, E. 1951. Theoretical Investigation Of An Automatic Control

System With Primary Sensitivity To Normal Accelerations As Used To Control A
Supersonic Canard Missile Configuration.

Seiler, P., Packard, A. and Gahinet, P. 2020. An Introduction to Disk Margins.

Tempo, R. and Dabbene, F. 2004. An Overview on Randomized Algorithms for Analysis

and Control of Uncertain Systems. Available from:
[Link]

Ujjin, R., Chaikiandee, S. and Ngaongam, C. 2021. Low Altitude Local Rocket
Aerodynamics Analysis and Experimental Testing In: basel, Switzerland: MDPI.

UKRA 2021. UK Rocketry Altitude Records | UKRA - United Kingdom Rocketry

Association. [Accessed 3 November 2022]. Available from:
[Link]

Uselton, J. 1971. Results of an Air Force Advanced Tactical Rocket Development

Program [Online]. [Accessed 1 May 2023]. Available from:
[Link]

Vigneswaran, R. and Kp, S. 2019. An Insight into the Dynamics and State Space
Modelling of a 3-D Quadrotor.

Yang, Y. 2012. Analytic LQR Design for Spacecraft Control System Based on
Quaternion Model. Journal of Aerospace Engineering. 25(3), pp.448–453.

Youds, T. 2022. Feasibility Analysis of a Reaction Wheel System for Sounding Rocket
Roll Control.

17
Youds, T., Cradock, B., Devlin, L., Whittaker, J. and Daney de Marcillac, J. 2023.
Development of a Canard Controlled Active Stability System for a Sounding
Rocket. Leeds, UK: University of Leeds.

18
 Appendix A – Supporting Formulae Derivation

This appendix has been lifted from Youds (2022) to support the derivation of the
aerodynamic equations in Chapter 3.

Normal Force Coefficient

The forces acting on the rocket in the body frame of reference can be divided into those
acting axially, i.e. along the 𝑥𝐵 axis, and those acting radially, i.e. normal to the 𝑥𝐵 axis.

The axial force components generally act through the rocket’s CG, or at such a small
offset that the resulting moments are much smaller than those generated by the normal
force components. Only the rocket’s rotational motion is of interest here and therefore
the axial force components have been omitted from the model.

Assuming that the thrust from the motor acts directly along the 𝑥𝐵 axis, the only normal
force components acting on the rocket are from gravity, when flying at a non-zero
pitch/yaw angle, and from aerodynamic lift and drag. Rockets fly at very small pitch and
yaw angles, and therefore the gravity component can generally be neglected with little
effect on the model’s accuracy.

Niskanen (2009) proposes that this aerodynamic normal force, 𝑁, can be calculated by
(A.1)
1
𝑁 = 𝜌𝑣 2 𝐴𝑟𝑒𝑓 𝐶𝑁 (𝐴. 1)
2
Where 𝜌 is the air density, 𝑣 is the freestream air velocity, 𝐴𝑟𝑒𝑓 is a reference area,
generally taken as the cross-sectional area of the rocket’s fuselage, and 𝐶𝑁 is the
rocket’s normal force coefficient.

This normal force is generated by the culmination of lift and drag generated when the
freestream air flows around the rocket body. Generally, rockets have axisymmetric body
components and fins with symmetrical aerofoils, and therefore generate no normal force
when flying at zero angle of attack, 𝛼 = 0. However, for 𝛼 > 0, these forces can be
significant and greatly affect the rotational motion of the rocket.

The variation of the normal force with angle of attack is described by 𝐶𝑁𝛼 , the angle of
attack derivative of the normal force coefficient – referred to as the normal force
coefficient derivative from here on. Niskanen (2009), extending on the work of
Barrowman (1967), states that the normal force coefficient varies linearly with angle of
attack, for small values of 𝛼. This relationship is given in equation (A.2).

19
 𝜕𝐶𝑁
𝐶𝑁𝛼 = | ⇒ 𝐶𝑁 ≈ 𝐶𝑁𝛼 𝛼 (𝐴. 2)
𝜕𝛼 𝛼=0
Niskanen (2009) provides methods for calculating the normal force coefficient derivative,
𝐶𝑁𝛼 , based on the rocket’s geometry and flight conditions, such as speed and angle of
attack. 𝐶𝑁𝛼 is calculated for each section of the rocket: nose, body and fins, and then
summed to find the value for the full vehicle.

Calculating 𝑪𝑵𝜶 for Axisymmetric Body Components

The nosecone and fuselage can be treated as two separate axisymmetric body
components. Niskanen (2009) derives equation (A.3) for finding 𝐶𝑁𝛼 of such
components. This approximation is very similar to that derived by Barrowman (1967) but
applies to a larger range of 𝛼 by not assuming that sin 𝛼 = 𝛼. The second term in
equation (A.3) is used to include the effect of body lift when flying at significant angles
of attack.

2 sin 𝛼 𝐴𝑝𝑙𝑎𝑛 sin2 𝛼

𝐶𝑁𝛼 = [𝐴(𝑙) − 𝐴(0)] +𝐾 (𝐴. 3)
𝐴𝑟𝑒𝑓 𝛼 𝐴𝑟𝑒𝑓 𝛼

Where 𝐴(𝑥) is the area at position 𝑥 along the component’s length, 𝐴𝑝𝑙𝑎𝑛 is the planform
area of the body component, and 𝐾 ≈ 1.1 is a correction factor. Care must be taken
when using this equation to avoid singularities at 𝛼 = 0. This is achieved by
implementing limit conditions, as shown by equations (A.4).

sin 𝛼 sin2 𝛼
lim =1 , lim =0 (𝐴. 4)
𝛼→0 𝛼 𝛼→0 𝛼
Calculating 𝑪𝑵𝜶 for Fins

Making use of Diederich’s semi-empirical method, Barrowman (1967) derived equation

(A.5) to calculate the normal force coefficient derivative for a single fin (Niskanen, 2009).
𝐴𝑓𝑖𝑛
𝐶 𝑁𝛼 𝐹𝐷 ( ) cos Γ𝑐
0 𝐴𝑟𝑒𝑓
(𝐶𝑁𝛼 ) = 4
(𝐴. 5)
1 2+𝐹𝐷 √1+ 2
𝐹𝐷

Where 𝐶𝑁𝛼0 is the normal force coefficient derivative of the 2D aerofoil, 𝐹𝐷 is Diedrich’s

planform correlation parameter, 𝐴𝑓𝑖𝑛 is the planform area of one fin, and Γ𝑐 is the
midchord sweep angle.

The 2D aerofoil normal force coefficient derivative, 𝐶𝑁𝛼0 , can be found based on thin

aerofoil theory, as in equation (A.6) (Niskanen, 2009).

20
 2𝜋
𝐶𝑁𝛼0 = (𝐴. 6)
√1−𝑀2

Where 𝑀 is the local Mach number. Gryphon I-C remains subsonic during flight,
therefore 𝑀 < 1.

Diederich provides equation (A.7) for calculating 𝐹𝐷 .

𝐴𝑅
𝐹𝐷 = 1 (𝐴. 7)
𝐶 cos Γ𝑐
2𝜋 𝑁𝛼0

Where 𝐴𝑅 is the fin aspect ratio. 𝐴𝑅 = 2𝑠 2 /𝐴𝑓𝑖𝑛 , where 𝑠 is the fin semi-span, can be
substituted into equation (A.7). Then substituting (A.7) and (A.6) into (A.5) gives an
equation for the normal force coefficient derivative for a single fin as a function of the
geometry and the Mach number. This is presented as equation (A.8).

𝑠2
2𝜋 𝐴
𝑟𝑒𝑓
(𝐶𝑁𝛼 ) = (𝐴. 8)
1 2
𝑠 2 √1 − 𝑀2
1 + √1 + ( 𝐴 cos Γ𝑐 )
𝑓𝑖𝑛

Niskanen (2009) then shows that the normal force coefficient derivative for 𝑘 fins is found
by multiplying the value for one fin by the sum of sin2 𝛬𝑖 where Λ𝑖 is the dihedral angle
of the fin relative to the incoming lateral airflow.
𝑘 𝑘

∑(𝐶𝑁𝛼 ) = (𝐶𝑁𝛼 ) ∑ sin2 Λ𝑖 (𝐴. 9)

𝑘 1
𝑖=1 𝑖=1

Fortunately, for rockets with three or more equally spaced fins, the sum term reduces to
𝑘/2. Thus, for Pathfinder, with four fins, equation (A.10) applies.

(𝐶𝑁𝛼 ) = 2(𝐶𝑁𝛼 ) (𝐴. 10)

4 1

Finally, a correction term, 𝐾 𝑇(𝐵) is applied to account for fin-body interference, defined
by equation (A.11).
𝑟𝑡
𝐾 𝑇(𝐵) = 1 + (𝐴. 11)
𝑠 + 𝑟𝑡

Where 𝑟𝑡 is the body radius at the fin position. Therefore, the final normal force coefficient
derivative for Pathfinder’s four fins is given by equation (A.12), where (𝐶𝑁𝛼 )1 is given by

equation (3.9).
𝑟𝑡
𝐶𝑁𝛼 𝑓𝑖𝑛𝑠 = 2(𝐶𝑁𝛼 )1 (1 + ) (𝐴. 12)
𝑠 + 𝑟𝑡

21
CP for Body Components

The centre of pressure of axially symmetric body components can be readily computed
by equation (A.13) (Niskanen, 2009).

𝑙 𝐴(𝑙 ) − 𝑉
𝑋𝐶𝑃 = (𝐴. 13)
𝐴(𝑙 ) − 𝐴(0)

Where 𝑉 is the volume of the component.

CP for Fins

The fin centre of pressure location can be obtained from an empirical equation originally
derived by Barrowman (1967).

𝑋𝑡 𝐶𝑟 + 2𝐶𝑡 1 𝐶𝑟2 + 𝐶𝑡2 + 𝐶𝑟 𝐶𝑡

𝑋𝐶𝑃 = + (𝐴. 14)
3 𝐶𝑟 + 𝐶𝑡 6 𝐶𝑟 + 𝐶𝑡

Where 𝑋𝑡 = 𝐶𝑟 − 𝐶𝑡 . This is measured from the leading edge so must be corrected to be

measured from the nose cone tip before being incorporated into finding the full rocket
CP.

Full Rocket CP

The full rocket centre of pressure is then found from the CP locations of each of the
individual components along with their respective 𝐶𝑁𝛼 values, as in equation (A.15)

∑𝑛𝑖=1 𝑋𝐶𝑃𝑖 (𝐶𝑁𝛼 )

𝑖
𝑋𝐶𝑃 = (𝐴. 15)
∑𝑛𝑖=1(𝐶𝑁𝛼 )
𝑖

22
 Appendix B – Linearisation

This appendix contains the mathematical derivation of the linearised dynamic equations
referenced in Chapter 4.

Roll

The roll dynamic equation is given by (B.1).

𝐽𝑅 𝑝̇ = 𝐿𝛼 [(𝛼 − 𝑥1 ) cos 𝑥1 − (𝛼 + 𝑥2 ) cos 𝑥2 + (𝛽 − 𝑦1 ) cos 𝑦1 − (𝛽 + 𝑦2 ) cos 𝑦2 ] − 𝐶𝑑𝑝 𝑝 (𝐵. 1)

Expanding the brackets gives (B.2).

𝐽𝑅 𝑝̇ = +𝐿𝛼 𝛼 cos 𝑥1
−𝐿𝛼 𝑥1 cos 𝑥1
−𝐿𝛼 𝛼 cos 𝑥2
−𝐿𝛼 𝑥2 cos 𝑥2
+𝐿𝛼 𝛽 cos 𝑦1 (𝐵. 2)
−𝐿𝛼 𝑦1 cos 𝑦1
−𝐿𝛼 𝛽 cos 𝑦2
−𝐿𝛼 𝑦2 cos 𝑦2
−𝐶𝑑𝑝 𝑝

The canard deflections are limited at ±10° and therefore it can be assumed that cos 𝑥1 =
cos 𝑥2 = cos 𝑦1 = cos 𝑦2 = 1, with only 1.5% error. Applying this assumption to (B.2),
collecting like terms and dividing through by 𝐽𝑅 gives (B.3), the linearised roll equation.

𝐿𝛼 𝐿𝛼 𝐿𝛼 𝐿𝛼 𝐶𝑑𝑝
𝑝̇ = − 𝑥1 − 𝑥2 − 𝑦1 − 𝑦2 − 𝑝 (𝐵. 3)
𝐽𝑅 𝐽𝑅 𝐽𝑅 𝐽𝑅 𝐽𝑅

Pitch

The pitch dynamic equation is given by (B.4).

𝐽𝐿 𝑞̇ + (𝐽𝑅 − 𝐽𝐿 )𝑝𝑟 = 𝑀𝛼 [(𝛼 − 𝑥1 ) cos 𝑥1 + (𝛼 + 𝑥2 ) cos 𝑥2 ] − 𝐶1 𝛼 − 𝐶2 𝑞 (𝐵. 4)

Rearranging and expanding the brackets gives (B.5).

23
 𝐽𝑅 𝑞̇ = +𝑀𝛼 𝛼 cos 𝑥1
−𝑀𝛼 𝑥1 cos 𝑥1
+𝑀𝛼 𝛼 cos 𝑥2
+𝑀𝛼 𝑥2 cos 𝑥2 (𝐵. 5)
−𝐶1 𝛼
−𝐶2 𝑞
−(𝐽𝑅 − 𝐽𝐿 )𝑝𝑟

The final term is not yet linear as it contains the product of two variables. Taylor series
linearisation can be applied to the final term, as in (B.6).

𝑓 (𝑝, 𝑟) = −(𝐽𝑅 − 𝐽𝐿 )𝑝𝑟

𝑑𝑓 𝑑𝑓
⇒ 𝑓 (𝑝, 𝑟) ≈ 𝑓(𝑝̅ , 𝑟̅ ) + |𝑝=𝑝̅ (𝑝 − 𝑝̅ ) + | (𝑟 − 𝑟̅ ) (𝐵. 6)
𝑑𝑝 𝑑𝑟 𝑟=𝑟̅

𝑝̅ and 𝑟̅ are steady state values, such that 𝑝̇ = 0 when evaluated at 𝑝 = 𝑝̅ and 𝑟 = 𝑟̅ . It
therefore follows that 𝑓 (𝑝̅ , 𝑟̅ ) = 0. Applying this simplification and substituting in deviation
variables, defined by (B.7), forms equation (B.8).

𝑝′ = 𝑝 − 𝑝̅ (𝐵. 7)
𝑟 ′ = 𝑟 − 𝑟̅
⇒ 𝑓(𝑝, 𝑟) ≈ −(𝐽𝑅 − 𝐽𝐿 )𝑟̅ 𝑝′ − (𝐽𝑅 − 𝐽𝐿 )𝑝̅ 𝑟 ′ (𝐵. 8)

There is an equilibrium in the system when all states in the state vector, 𝐱, equal zero.
Therefore 𝑝̅ = 𝑟̅ = 0 and consequently 𝑓(𝑝, 𝑟) = 0.

Substituting this back into (B.6), applying the small angle approximation to the cosines
in the first four terms, and dividing through by 𝐽𝐿 gives (B.9), the linearised pitch motion
equation.
𝑀𝛼 𝑀𝛼 𝑀𝛼 𝐶1 𝐶2
𝑞̇ = − 𝑥1 + 𝑥2 + (2 − )𝛼 − 𝑞 (𝐵. 9)
𝐽𝐿 𝐽𝐿 𝐽𝐿 𝐽𝐿 𝐽𝐿

Yaw

The nonlinear yaw dynamic equation is given by (B.10).

𝐽𝐿 𝑟̇ + (𝐽𝐿 − 𝐽𝑅 )𝑝𝑞 = 𝑁𝛼 [(𝛽 − 𝑦1 ) cos 𝑦1 + (𝛽 + 𝑦2 ) cos 𝑦2 ] − 𝐶1 𝛽 − 𝐶2 𝑟 (𝐵. 10)

The same process used to linearise the pitch equation was used for the yaw equation in
(B.10). This gives (B.11), the linearised yaw motion equation.
𝑁𝛼 𝑁𝛼 𝑁𝛼 𝐶1 𝐶2
𝑟̇ = − 𝑦1 + 𝑦2 + (2 − )𝛽 − 𝑟 (𝐵. 11)
𝐽𝐿 𝐽𝐿 𝐽𝐿 𝐽𝐿 𝐽𝐿

24
State Space Formulation

Equations (B.3), (B.9), and (B.11) provide a linear approximation of the rotational motion
of the rocket in terms of the states, 𝐱, inputs, 𝐮, and disturbances, 𝐝. This linear system
of equations can be written as matrices in state space form, as in (B.12).
0 0 0 1 0 0 0 0 0 0
0 0 0 0 1 0 0 0 0 0 0 0
𝑝 0 0 0 0 0 1 𝜙 0 0 0 0 0 0
𝑞 𝐶𝑑𝑝 𝜃 𝐿𝛼 𝐿𝛼 𝐿𝛼 𝐿𝛼 𝑥1 0 0
𝑟 0 0 0 − 0 0 𝜓 − − − − 𝑥 0 0 𝛼
𝐽𝑅 𝐽𝑅 𝐽𝑅 𝐽𝑅 𝐽𝑅 { 2 } + 𝑀 {𝛽} (𝐵. 12)
𝑝̇ = 𝐶2 𝑝
+
𝑀𝛼 𝑀𝛼 𝑦1
2
𝛼
−
𝐶1
0
𝑞̇ 0 0 0 0 − 0 𝑞 − + 0 0 𝑦2 𝐽𝐿 𝐽𝐿
{𝑟̇ } 𝐽𝐿 {𝑟} 𝐽𝐿 𝐽𝐿 𝑁𝛼 𝐶1
𝐶2 𝑁𝛼 𝑁𝛼 0 2 −
0 0 0 0 0 − 0 0 − + [ 𝐽𝐿 𝐽𝐿 ]
[ 𝐽𝐿 ] [ 𝐽𝐿 𝐽𝐿 ]

25
 Appendix C – Expected Parameter Variations

Table 2 displays the expected variations of parameters within the nonlinear model.
These were generally obtained from taking the worst-case tolerance stack-up from the
component’s respective manufacturing processes. A very conservative approach was
taken for aerodynamic components, especially those manufactured by hand, such as
the rear fins.

The ±50% variation on the canard lift coefficient derivative is to account for the
inaccuracy of the aerodynamic simulations.
Table 2 - Parameter variations used for the robustness analysis of the control system.

Parameter Nominal Value Variation Units

Radial moment of inertia 𝐽𝑅 * ± 20% kg m2
Longitudinal moment of inertia 𝐽𝐿 * ± 20% kg m2
Rocket body diameter 𝑑𝑏𝑜𝑑𝑦 0.103 ± 0.003 m
Nose cone length 𝑙𝑛𝑜𝑠𝑒 0.507 ± 0.010 m
Rocket body length 𝑙𝑏𝑜𝑑𝑦 1.563 ± 0.010 m
Fin semi-span 𝑠𝑠𝑓𝑖𝑛 0.14 ± 0.002 m
Fin root chord length 𝑐𝑟𝑓𝑖𝑛 0.12 ± 0.002 m
Fin tip chord length 𝑐𝑡𝑓𝑖𝑛 0.065 ± 0.002 m
Rocket CG location 𝑋𝐶𝐺 1.32 ± 0.050 m
Canard CP location (𝑥) 𝑋𝐶𝑃𝑐𝑎𝑛 0.589 ± 0.010 m
Canard CP location (𝑦) 𝑌𝐶𝑃𝑐𝑎𝑛 0.0865 ± 0.005 m
Canard CP location (𝑧) 𝑍𝐶𝑃𝑐𝑎𝑛 0.0865 ± 0.005 m
Canard semi-span 𝑠𝑠𝑐𝑎𝑛 0.07 ± 0.001 m
Canard root chord length 𝑐𝑟𝑐𝑎𝑛 0.06 ± 0.001 m
Canard tip chord length 𝑐𝑡𝑐𝑎𝑛 0.05 ± 0.001 m
Canard lift coefficient derivative 𝐶𝐿𝛼 2.86 ± 50% rad-1

* 𝐽𝑅 and 𝐽𝐿 vary during flight, and their profiles were exported from OpenRocket. The
variation of ±20% was applied to the full profile before importing into Simulink. This
maintained the same variation during flight, but shifted it either up or down, depending
on the value used for that simulation.

26
 Appendix D – Canard Moment Diagrams

Figure 5 shows the geometry of the pitch and yaw canard configuration with respect to
the CG, to aid with understanding the derivation of the dynamics equations in Chapter
3.

This diagram shows the 𝑥1 canard. The 𝑥2 canard is hidden behind the rocket body.

Figure 5 - Pitch / yaw canard configuration.

Figure 6 shows the configuration for the roll axis.

Figure 6 - Roll canard configuration.

27
 Appendix E – Enhanced Wind Modelling Results

Figure 7 presents the system response to a varying crosswind both with and without the
control system enabled. Dryden’s continuous wind turbulence model was implemented
in Simulink (MathWorks, n.d.), where both the speed and direction of the wind vector
varies with the rocket’s altitude and attitude.

Figure 7 - System performance in response to a varying wind vector.

28