Open access peer-reviewed chapter

# State-Space Modeling of a Rocket for Optimal Control System Design

Written By

Aliyu Bhar Kisabo and Aliyu Funmilayo Adebimpe

Reviewed: October 29th, 2018 Published: June 5th, 2019

DOI: 10.5772/intechopen.82292

From the Edited Volume

## Ballistics

Edited by Charles Osheku

Chapter metrics overview

2,439 Chapter Downloads

View Full Metrics

## Abstract

This chapter is the first of two others that will follow (a three-chapter series). Here we present the derivation of the mathematical model for a rocket’s autopilots in state space. The basic equations defining the airframe dynamics of a typical six degrees of freedom (6DoFs) are nonlinear and coupled. Separation of these nonlinear coupled dynamics is presented in this chapter to isolate the lateral dynamics from the longitudinal dynamics. Also, the need to determine aerodynamic coefficients and their derivative components is brought to light here. This is the crux of the equation. Methods of obtaining such coefficients and their derivatives in a sequential form are also put forward. After the aerodynamic coefficients and their derivatives are obtained, the next step is to trim and linearize the decoupled nonlinear 6DoFs. In a novel way, we presented the linearization of the decoupled 6DoF equations in a generalized form. This should provide a lucid and easy way to implement trim and linearization in a computer program. The longitudinal model of a rocket presented in this chapter will serve as the main mathematical model in two other chapters that follow in this book.

### Keywords

• rocket
• six degrees of freedom (6DoFs)
• state space
• trimming
• linearization

## 1. Introduction

Over the years, a number of authors [1, 2, 3, 4] have considered modeling rocket/missile dynamics for atmospheric flights. In the majority of the published work on these mathematical models, trimming and locally linearization are done without detailed explanation to the variables in the decoupled airframe dynamics. As such, the easy-to-write computer programs to facilitate this process numerically have been impeded.

With the advent of fast processors and numerical software like MATLAB®, Maple, Python, etc., it is now possible to take a complex nonlinear 6DoF equation like that of a rocket and run a program that can trim and linearize it with ease. This has made the field of control system design to grow at an exponential rate [5].

It is a known fact that the mathematical models are developed with their use in mind. This means before delving into the realization of a model, one should be well informed of the purpose for which that mathematical model will serve [6]. Especially, in the field of control system design, a mathematical model in transfer function might not be ideal for optimal control design. However, problem-solving environments (PSEs) like MATLAB®/Simulink come with built-in functions capable of transforming, say state-space model, to transfer functions. One should bear in mind that this is not without a “cost.”

Advertisement

## 2. Mathematical model

A nontrivial part of any control problem is modeling the process. The objective is to obtain the simplest mathematical description that adequately predicts the response of the physical system to all inputs. For a rigid dynamic body, its motion can be described in translational, rotational, and angular inclinations at all times.

### 2.1. Translational motion

An accelerometer is often used to measure force on a dynamic body. For a rocket in motion, these forces [7] are represented as given in (1). Actually, this is a measure of the specific force, i.e., the nongravitational force per unit mass in x,y,z-directions, respectively. The specific force (also called the g-force or mass-specific force) has units of acceleration or ms−1. So, it is not actually a force at all but a type of acceleration:

u ̇ = F A x b + F P x b + F g x b m qw rv , m / s 2 v ̇ = F A y b + F P y b + F g y b m ru pw , m / s 2 w ̇ = F A z b + F P z b + F g z b m pv qu , m / s 2 E1

where

F A xb , F A yb , F A zb = components of aerodynamic force vector FA expressed in the body.

coordinate system, N.

F g xb , F g yb , F g zb = components of gravitational force vector Fg expressed in body coordinate system, N.

F p xb , F p yb , F p zb = components of thrust vector Fp expressed in the body coordinate system, N.

m = instantaneous rocket mass, kg.

p, q, r = components of angular rate vector ω expressed in body coordinate system (roll, pitch, and yaw, respectively), rad/s.

u, v, w = components of absolute linear velocity vector Vx expressed in the body coordinate system, m/s.

u ̇ , v ̇ , w ̇ = components of linear acceleration expressed in body coordinate system, m/s2.

The aerodynamic forces have the following components:

F A x b = 0.5 ρ V M 2 C A S F A y b = 0.5 ρ V M 2 C N y S F A z b = 0.5 ρ V M 2 C N z S E2

where

CA = aerodynamic axial force coefficient, dimensionless.

CN = aerodynamic normal force coefficient, dimensionless.

CNy = coefficient corresponding to component of normal force on yb-axis

C N y = C N v v 2 + w 2 , dimensionless

CNz = coefficient corresponding to component of normal force on zb-axis

C N z = C N 1 v 2 + w 2 , dimensionless

S = aerodynamic reference area, m2.

VM = magnitude of velocity vector of the center of mass of the rocket, m/s.

ρ = atmospheric density, kg/m3.

The propulsive forces in (1), as depicted in (Figure 1), are computed [8] as follows:

F p x b = F p cos γ 2 cos γ 1 , N F p y b = F p cos γ 2 sin γ 1 , N F p z b = F p sin γ 2 , N E3

where

Fp = magnitude of total instantaneous thrust force vector, N.

Fpxb, Fpyb, Fpzb = components of thrust vector Fp expressed in body coordinate system, N.

γ1 = angle measured from xb-axis to projection of thrust vector Fp on xb yb-plane, rad (deg).

γ2 = angle measured from projection of thrust vector Fp on xb yb-plane to the thrust vector Fp, rad.

lp = the distance from the aerodynamic center to center of mass, m.

l = the distance from the center of mass to nozzle, m.

If the rocket was designed for thrust vector control via gimbaling of the nozzle, Fp will be computed as given in (3). Here we assume that Fp is acting along the line of symmetry of the rocket because the nozzle is fixed (fin control). Hence, the magnitude of the thrust force Fp is calculated by

F p = F pref + p ref p a A e , N E4

where

Ae = rocket nozzle exit area, m2.

Fpref = reference thrust force, N.

Pa = ambient atmospheric pressure, Pa.

Pref = reference ambient pressure, Pa.

The gravitational forces in (1) are computed as follows:

F g x e = 0 , N F g y e = 0 , N F g z e = mg , N E5

where

Fgxe, Fgye, Fgze = components of gravitational force vector Fg expressed in earth coordinate system, N.

g = acceleration due to gravity, m/s2.

m = instantaneous mass of rocket, kg.

The dependence of the acceleration due to gravity on the altitude of the rocket is given by

g = g 0 R e 2 R e + h 2 , m / s 2 E6

where

g = acceleration due to gravity, m/s2.

g0 = acceleration due to gravity at earth surface (nominally 9.8 m/s2), m/s2.

h = altitude above sea level, m.

Re = radius of the earth, m.

The gravitational force expressed in body coordinates is computed by multiplying (5) by matrix (7):

T e / b = cθcψ sϕsθcψ cϕsψ cϕsθcψ + sϕsψ cθsψ sϕsθsψ + cϕcψ cϕsθcψ sϕcψ sϕcθ cϕcθ , dimensionless E7

where

c = cosine function (cθ = cos θ), dimensionless.

s = sine function (sθ = sin θ), dimensionless.

θ = the Euler angle of rotation in elevation (pitch) of body frame relative to earth frame, rad (deg).

φ = the Euler angle of rotation in roll of body frame relative to earth frame, rad (deg).

ψ = the Euler angle of rotation in azimuth (heading) of body frame relative to earth frame, rad (deg).

[Te/b] = transformation matrix from earth to body.

A vector v expressed in the body coordinate system can be transformed to the earth coordinate system by the matrix equation

v e = T e / b v b E8

Hence, considering (5) we can write the following:

F g x b F g y b F g z b = T b / e F g x e F g y e F g z e , N E9

The terms Fgxb, Fgyb, and Fgzb are the components of the gravitational force substituted into (1).

The mass in (1) is given below:

m = m 0 1 I sp 0 t F p ref dt , kg E10

where

Fpref = reference thrust force, N.

Isp = specific impulse of propellant, Ns/kg.

m0 = rocket mass at time zero (i.e., at the time of launch), kg.

t = simulated time, s.

### 2.2. Rotational motion

A gyroscope or gyro is a device that measures the angular acceleration or rotational motion of a dynamic body. On a rocket, this rotational motion can be described as

p ̇ = L A + L p qr I z I y I x , rad / s 2 deg / s 2 q ̇ = M A + M p rp I x I z I y , rad / s 2 deg / s 2 r ̇ = N A + N p pq I y I x I z , rad / s 2 deg / s 2 E11

where

LA, MA, NA = components of aerodynamic moment vector MA expressed in body coordinate system (roll, pitch, and yaw, respectively), Nm.

Lp, Mp, Np = components of propulsion moment vector Mp expressed in body coordinate system (roll, pitch, and yaw, respectively), Nm.

Ix, Iy, Iz = components of inertia (diagonal elements of inertia matrix when products of inertia are zero), kg-m2.

p, q, r (P, Q, R= components of angular rate vector ω expressed in body coordinate system (roll, pitch, and yaw, respectively), rad/s (deg/s).

p ̇ , q ̇ , r ̇ = components of angular acceleration ω ̇ expressed in body coordinate (roll, pitch, and yaw, respectively), rad/s2 (deg/s2).

L A = 0.5 ρ V M 2 C l Sd , Nm M A = 0.5 ρ V M 2 C m Sd , Nm N A = 0.5 ρ V M 2 C n Sd , Nm E12

where

Cl = aerodynamic roll moment coefficient about center of mass, dimensionless.

Cm = aerodynamic pitch moment coefficient about center of mass, dimensionless.

Cn = aerodynamic yaw moment coefficient about center of mass, dimensionless.

d = aerodynamic reference length of body, m.

The aerodynamic moment coefficients are of the form

C l = C l δ δ r + d 2 V M C l p p C m = C m ref C N z x cm x ref d + d 2 V M C m q + C m α q C n = C n ref + C N y x cm x ref d + d 2 V M C n r + C n β r dimensionless E13

where

Clp = roll damping derivative relative to roll rate p, rad−1 (deg−1).

C = slope of curve formed by roll moment coefficient C1 versus control surface deflection, rad−1(deg−1).

Cmref = pitching moment coefficient about reference moment station, dimensionless.

Cmq = pitch damping derivatives relative to pitch rate q, rad−1(deg−1).

C m α ̇ = pitch damping derivative relative to angle of attack rate α ̇ (slope of curve formed by Cα verses α ), rad−1 (deg−1).

CNy = coefficient corresponding to component of normal force on yb-axis, dimensionless.

CNz = coefficient corresponding to component of normal force on zb-axis, dimensionless.

C n r ̇ = yaw damping derivative relative to yaw rate r ̇ , rad−1(deg−1).

Cnref = yawing moment coefficient about reference moment station, dimensionless.

C n β ̇ = yaw damping derivative relative to angle of sideslip rate β ̇ , rad-l (deg-1).

xcm = instantaneous distance from rocket nose to center of mass, m.

xref = distance from rocket nose to reference moment station, m.

δr = effective control surface deflection causing rolling moment, rad (deg).

Considering that the rocket we intend to control is via fin deflection, fins on the rocket will be designated as shown in Figure 2.

Hence, the following moment coefficients are also given as

C m ref = C m α α + C m δ δ η C n ref = C n β β + C n δ δ ζ , dimensionless E14

where

Cmref = pitching moment coefficient about reference moment station (this is the static value normally measured in the wind tunnel.), dimensionless.

C m α = slope of curve formed by pitch moment coefficient. Cm versus angle of attack, α rad−1 (deg−1) slope of tune formed.

C m δ = slope of tune formed by pitch moment coefficient Cm versus control surface deflection for pitch, δp rad−1 (deg−1).

C n β = slope of curve formed by yawing moment coefficient Cn versus angle of sideslip β, rad−1 (deg−1).

C n δ = slope of curve formed by yaw moment coefficient Cn versus effective control surface deflection for yaw, δy rad−1 (deg−1).

α = angle of attack, rad (deg).

β = angle of sideslip, rad (deg).

δη = δ1 = δ3 effective control surface deflection causing pitching moment, rad (deg).

δζ = δ4 = δ2 = effective control surface deflection causing yawing moment, rad (deg).

The angle of attack, angle of sideslip, and roll angle required for the realization of the aerodynamic coefficients are

α = Tan 1 w u or α t = tan 1 v 2 + w 2 u β = sin 1 v V M ϕ = tan 1 v w , rad deg E15

where ∅ is aerodynamic roll angle, rad (deg), and αt is the total angle of attack.

Table 1 gives a list of the aerodynamic coefficients that must be obtained for every rocket design before a model can be realized. There exist numerical and semi-numerical means of obtaining such coefficients. Examples of software that can do semiempirical computation of such coefficients and their derivatives are Missile Digital DATCOM® [9] and Flexible Structures Simulator (FSS) [10]. Finally, a wind tunnel test is expected to validate and update all coefficients and their derivatives before the system engineer delves in the control design.

S/N Variable Description
1 CN Normal force coefficient (body axis)
2 CL Lift coefficient (wind axis)
3 CM Pitching moment coefficient (body axis)
4 Xcp Center of pressure in calibers from moment reference center
5 CA Axial force coefficient (body axis)
6 CD Drag coefficient (wind axis)
7 CY Side force coefficient (body axis)
8 Cn Yawing moment coefficient (body axis)
9 Cl Rolling moment coefficient (body axis)
10 C Normal force coefficient derivative with angle of attack
11 C Pitching moment coefficient derivative with angle of attack
12 C Side force coefficient derivative with sideslip angle
13 C Yawing moment coefficient derivative with sideslip angle (body axis)
14 C Rolling moment coefficient derivative with sideslip angle (body axis)
15 CMq Pitching moment coefficient derivative with pitch rate
16 CNq Normal force coefficient derivative with pitch rate
17 CAq Axial force coefficient derivative with pitch rate
18 C M α ̇ Pitching moment derivative with rate of change of angle of attack
19 C N α ̇ Normal force derivative with rate of change of angle of attack
20 Clp Rolling moment coefficient derivative with roll rate
21 Cnp Yawing moment coefficient derivative with roll rate
22 CYp Side force coefficient derivative with roll rate
23 Clr Rolling moment coefficient derivative with yaw rate
24 Cnr Yawing moment coefficient derivative with yaw rate
25 CYr Side force coefficient derivative with yaw rate

### Table 1.

Aerodynamic coefficients as a function of angle of attack.

The third and final component to fully describe the motion of a rocket is its angular inclination or attitude. We chose the Euler angles to describe the attitude of the rocket.

### 2.3. The Euler angles

Missile attitude is required for a number of simulation functions including the calculation of angle of attack, seeker gimbal angles, fuze look angles, and warhead spray pattern. In simulations with six degrees of freedom, the missile attitude is calculated directly by integrating the set of equations that define the Euler angle rates, i.e.,

ϕ ̇ = p + q sin ϕ + r cos ϕ tan θ , rad / s deg / s θ ̇ = q cos ϕ r sin ϕ , rad / s deg / s ψ ̇ = q sin ϕ + r cos ϕ cos θ , rad / s deg / s E16

where

θ = the Euler angle rotation in elevation (pitch angle), rad (deg).

ϕ = the Euler angle rotation in roll (roll angle), rad (deg).

ψ ̇ , θ ̇ , ϕ ̇ = rates of change of the Euler angles in heading, pitch, and roll, respectively, rad/s (deg/s).

P, q, r = components of angular rate vector w expressed in body coordinate system (roll, pitch, and yaw, respectively), rad/s.

The three heading angle of ψ ̇ , θ ̇ , ϕ ̇ can be measured directly using a magnetometer. Note also that the Euler angles in (16) can all be derived by integrating gyroscopic measurements. As such we might not need an instrument that will measure it directly. Nevertheless, a magnetometer can be added to the instrumentation on board to measure heading.

Combining (1), (11), and (16) gives a complete six-degree-of-freedom equation of motion for a rocket in flight as shown in Figure 3. This could be written together as given in (17). In today’s modern aerospace industry, a single device like the MPU6050, a MEM-based integrated chip, can be used to give numerical values for state variables of (17) on any dynamic body. For control system design, the rocket system as described in (17) needs to be separated into the two planes (decouple); these are called the lateral (la) and longitudinal (lo) dynamic equations of motion:

u ̇ = F A x b + F P x b + F g x b m qw rv , m / s 2 v ̇ = F A y b + F P y b + F g y b m ru pw , m / s 2 w ̇ = F A z b + F P z b + F g z b m pv qu , m / s 2
p ̇ = L A + L p qr I z I y I x , rad / s 2 deg / s 2 q ̇ = M A + M p rp I x I z I y , rad / s 2 deg / s 2 r ̇ = N A + N p pq I y I x I z , rad / s 2 deg / s 2 E17
ϕ ̇ = p + q sin ϕ + r cos ϕ tan θ , rad / s deg / s θ ̇ = q cos ϕ r sin ϕ , rad / s deg / s ψ ̇ = q sin ϕ + r cos ϕ cos θ , rad / s deg / s

Since the control system we intend to design is in a family of linear controllers and (17) is a nonlinear system of differential equations, trimming and linearization must be done.

### 2.4. Trimming

To trim or find equilibrium values requires a good knowledge of advance computational techniques. A trim point, also known as an equilibrium point, is a point in the parameter space of a dynamic system at which the system is in a steady state. The trim problem for a rocket can be described as finding a set of suitable input values to satisfy a set of conditions. Hence, a trim point involves setting of its controls that causes the rocket to fly straight and level in the longitudinal plane. The suitable input values are the control surface deflections, the thrust setting, and the rocket attitude [11]. The set of conditions are the rocket’s accelerations. The variables associated with the trim problems can be divided into three categories:

• Objective variables

• Control variables

• Flight condition variables

The objective variables need to be driven toward the specified values, often zero (i.e., steady flight with zero sideslip). The objective parameters are combined in the objective vectors ola (state vectors) and olo, for the lateral and longitudinal missile dynamics as

o la = v ̇ p ̇ r ̇ β T . E18
o lo = u ̇ w ̇ q ̇ α T . E19

The sideslip angle is also included, since for most cases, there are multiple solutions to the trim problem, each with a different sideslip angle. In the desired solution, the sideslip angle should be zero. In that case, the drag is at a minimum. The control parameters are adjusted in order to drive the objective parameters to their specified values. Together, they form the control vectors cla and clo, described in (22). The control variables (input variables) describe the trimmed pilot control and the aircraft attitude:

c la = δ a δ r ϕ ψ T , E20
c lo = δ e τ θ T , E21

Finally, the 12 states of the 6DoF equation of motion must be initialized with the initial state conditions. In MATLAB®, the trim command is used to find equilibrium points. The object of trimming is to bring the forces and moments acting on the rocket into a state of equilibrium. That is the condition when the axial, normal, and side forces and the roll, pitch, and yaw moments are all zero. Mathematically, trimming combines implicit and explicit Jacobian approaches. The Jacobian trim approach is the preferred method for rigid rocket. The Jacobian method is robust, since trim convergence is likely to occur even with rough estimates of the Jacobian and a rough first guess. Note that in general, any optimization routine could be used to solve the trim problem, as long as it is robust enough.

The Jacobian approach is based on the assumption that the change of the objective vector is linearly related to the change in the control input, which is shown in (22):

o i + 1 o i = J i c i + 1 c i . E22

In (22), Ji is the Jacobian matrix evaluated near control input ci. Its entries are first-order partial derivatives and represent the effect of changes in each control input on acceleration.

Note that changes in lateral plane induce changes in the longitudinal plane and vice versa; thus, we can write (23) the Jacobian for the lateral dynamics and (24) for the longitudinal or pitch dynamics:

J i = o c c i la = u ̇ δ a u ̇ δ r u ̇ ϕ u ̇ ψ v ̇ δ a v ̇ δ r v ̇ ϕ v ̇ ψ w ̇ δ a w ̇ δ r w ̇ ϕ w ̇ ψ p ̇ δ a p ̇ δ r p ̇ ϕ p ̇ ψ q ̇ δ a q ̇ δ r q ̇ ϕ q ̇ ψ r ̇ δ a r ̇ δ r r ̇ ϕ r ̇ ψ β δ a β δ r β ϕ β ψ c i . E23
J i = o c c i lo = u ̇ δ e u ̇ τ u ̇ θ v ̇ δ e v ̇ τ v ̇ θ w ̇ δ e w ̇ τ w ̇ θ p ̇ δ e p ̇ τ p ̇ θ q ̇ δ e q ̇ τ q ̇ θ r ̇ δ e r ̇ τ r ̇ θ β δ e β τ β θ c i . E24

If the rocket is assumed to be at equilibrium, or trim condition, then the equations of motion can be linearized, and the 6DoF equation of motion can be resolved into their lateral and longitudinal states.

### 2.5. Linearization

The system of first-order nonlinear differential equations of the rocket as presented by (17) is said to be in state variable form if its mathematical model is described by a system of n first-order differential equations and an algebraic output equation as [12]:

x ̇ 1 = f 1 x 1 x n u x ̇ 2 = f 2 x 1 x n u x ̇ n = f n x 1 x n u y = h x 1 x n u E25

The column vector x = [x1,…xn]T is called the state of the system. The scalars u and y are called the control input and the system output, respectively, denoting

f x u = f 1 x 1 x n u f 2 x 1 x n u f n x 1 x n u . E26

Concisely, (26) is written as

x ̇ = f x u , y = h x u . E27

where f and h are nonlinear functions of x and u; then, we say that the system is nonlinear. To linearize (26), we desire it to become

x ̇ = Ax + Bu y = Cx + Du . E28

where A is n x n, B is n x 1, C is 1 x n, and D is all scalar in MATLAB®/Simulink.

One reason for approximating the nonlinear system (26) by a linear model of the form (28) is that, by so doing, one can apply rather simple and systematic linear control design techniques. Given the nonlinear system (26) and an equilibrium or trimmed points x = [ x 1 x n ]T obtained when u = u, noting that x = x = x , we define a coordinate transformation as follows:

Δx = Δ x 1 Δ x n = x 1 x 1 x n x n .

Further, denoting u = u u , y = y h x u . The new coordinates x , u , and y represent the variations of x, u, and y from their equilibrium values. You have to think of these as a new state, new control input, and new output, respectively. The linearization of (26) is at x∗ in which the equilibrium or trim [12] state is given by

Δ x ̇ = AΔx + BΔu Δy = CΔx + DΔu , E29

where

A = f u x , u = f 1 x 1 x 1 x n u f 1 x n x 1 x n u f n x 1 x 1 x n u f n x n x 1 x n u ,
B = f u x , u = f 1 u x 1 x n u f n u x 1 x n u ,
C = h x x , u = h x 1 x 1 x n u h x n x 1 x n u , D = h u x , u .

Note that the matrices A, B, C, and D have constant coefficients in that all partial derivatives are evaluated at the numerical values x 1 x n u .

### 2.6. Lateral dynamics of a rocket

Equations of motion in the lateral plane are described by (30). Note that (30) comprises of one of the force equations (Fy), one of the momentum equations (My), and two of the Euler angles from the 6DoF equations (decoupled from (17)):

v ̇ = F A y b + F P y b + F g y b m ru pw , m / s 2 p ̇ = L A + L p qr I z I y I x , rad / s 2 deg / s 2 r ̇ = N A + N p pq I y I x I z , rad / s 2 deg / s 2 ϕ ̇ = p + q sin ϕ + r cos ϕ tan θ , rad / s deg / s ψ ̇ = q sin ϕ + r cos ϕ cos θ , rad / s deg / s E30

For a completely computed aerodynamic coefficients and their derivatives, (30) can be expressed in state-space form [13] as

v ̇ p ̇ r ̇ ϕ ̇ ψ ̇ = y v y p y r y ϕ y ψ l v l p l r l ϕ l ψ n v n p n r n ϕ n ψ 0 1 0 0 0 0 0 1 0 0 v p r ϕ ψ + y ξ y ζ l ξ l ζ n ξ n ζ 0 0 0 0 ξ ζ . E31

All the variables in A matrix of (31) are the lateral dimensionless aerodynamic stability derivatives with respect to the system state vectors. The variables in B matrix of (31) are the lateral dimensionless aerodynamic control derivatives with respect to the designated control surfaces.

### 2.7. Longitudinal dynamics of a rocket

The longitudinal dynamics in motion is also called the pitch plane. Equations describing the motion of a rocket in this plane can be described as given in (32). Note that (32) comprises two of the force equations (Fx and Fz), two of the momentum equations (Mx and Mz), and two of the Euler angles from the 6DoF equations as given in (17):

u ̇ = F A x b + F P x b + F g x b m qw rv , m / s 2 w ̇ = F A z b + F P z b + F g z b m pv qu , m / s 2 q ̇ = M A + M p rp I x I z I y , rad / s 2 deg / s 2 θ ̇ = q cos ϕ r sin ϕ , rad / s deg / s E32

Just as with (31), (32) can also be re-expressed in state space as

u ̇ w ̇ q ̇ θ ̇ = x u x w x q x θ z u z w z q z θ m u m w m q m θ 0 0 1 0 u w q θ + x η x τ z η z τ m η m τ 0 0 η τ E33

In MATLAB® the linmod [14] command is used to invoke linearization. The assumption made for decoupling the linear model is that the cross-coupling effects between the two modes are negligible. These assumptions are:

• The rocket is designed with conventional control surfaces that do not give significant cross-coupling control between the lateral and longitudinal modes.

• The rocket is symmetrical about the xz plane in which the inertia cross-coupling in (xy and xz planes) result to cross-coupling between the lateral and longitudinal modes is minimum.

It can be shown that a typical trimmed and linearized model of the pitch plane (longitudinal dynamics) for a rocket [15] is given as presented in (34). Notice that compared to (33), the velocity in x-direction is not included in (34). This is basically due to the fact that in this pitch plane, translational motion for a rocket is predominantly in the z-direction (velocity w):

θ ̇ q ̇ w ̇ = 0 1 0 14.7805 0 0.01958 100.858 1 0.1256 θ q w + 0 3.4858 20.42 δ η + 0 14.7805 94.8557 α w y = 1 0 0 θ q w E34

where

A = 0 1 0 14.7805 0 0.01958 100.858 1 0.1256 , B = 0 3.4858 20.42 , G = 0 14.7805 94.8557 C = 1 0 0 , D = 0 .
Advertisement

## 3. Discussion of result

From (34), it can be seen that a three state variable models have been realized in state space. This implies that modern observer like the Kalman filter can be designed to estimate and predict the trajectory of such rocket dynamics. This mathematical model also can be used to design all the control algorithms that fall in the class of modern (optimal theory) control. Particularly, this model is important in the realization of the longitudinal autopilot system of a rocket.

Advertisement

## 4. Conclusion

Mathematical models are the bedrock of almost all scientific activities. Here we were able to define the nonlinear airframe dynamics of a rocket in 6DoF. We went further to decouple the 6DoF equations of motion for the rocket and presented forms in which the decoupled 6DoF is easily trimmed and linearized with a computer program like MATLAB®. The process presented here can be repeated for any size of the rocket and aircrafts/unmanned aerial vehicle (UAV). Note that if the aerospace vehicle being model is not a rocket, and a state-space model is needed, all the procedures outlined in this chapter will be the same. The only changes that will be accommodated will come from the numerical values of the aerodynamic coefficients and their derivatives.

## References

1. 1. Papp I. Missile mathematical model and system design. Atlantic Association for Research in the Mathematical Sciences. 2017;16(1):29-35
2. 2. Farhan AF. State space model for autopilot design of aerospace vehicles. In: Weapons Systems Division. Defence Science and Technology Organisation, DSTO-TR-1990; 2007
3. 3. Chelaru T-V, Barbu C. Mathematical model for sounding rockets, using attitude and rotation angles. International Journal of Applied Mathematics and Informatics. 2009;3(2):35-44
4. 4. Belega B-A. A mathematical model for computing the trajectory of rockets in a resistant medium taking into account the Earth’s rotation. The system differential equations of the rocket trajectory. In: International Conference of Scientific Paper AFASES 2012; 2012
5. 5. Siouris GM. Missile Guidance and Control Systems. New York: Springer-Verlag, Inc; 2004. ISBN: 0-387-00726-1
6. 6. Blakelock JH. Automatic Control of Aircraft and Missiles. New York: John Wiley & Sons, Inc; 1965
7. 7. Missile Flight Simulation Part One Surface-to-Air Missiles, MIL-HDBK-1211(MI) 17 July 1995. In: Military Handbook. 1995
8. 8. Tewari A. Advance Control of Aircraft, Spacecraft and Rockets. United Kingdom: John Wiley & Sons, Ltd; 2011. ISBN 9781119971207
9. 9. MISSILE DATCOM User’s Manual—2011 Revision. Christopher Rosema, Joshua Doyle, Lamar Auman, and Mark Underwood-U.S. Army Aviation & Missile Research, Development and Engineering Centre. William B. Blake-Control Design and Analysis Branch Control Sciences Division. AFRL-RB-WP-TR-2011-3071 (Supersedes 2008 Revision, AFRL-RB-WP-TR-2009-3015)
10. 10. User Manual. Software Package for Flexible Structures Simulation (FSS). International Institute for Advanced Aerospace Technologies of St. Petersburg State University of Aerospace Instrumentation. Performed according to the Contract No. I-011-2012; 2013
11. 11. Yang WY et al. Applied Numerical Methods Using MATLAB. Hoboken, New Jersey: John Wiley & Sons, Inc; 2005. ISBN 0-471-69833-4
12. 12. White RE. Computational Mathematics: Models, Methods, and Analysis with MATLAB and MPI. New York: Chaoman & Hall/CRC; 2004. ISBN: 1-58488-364-2
13. 13. Cook MV. Flight Dynamics Principles: A Linear Approach to Aircraft Stability and Control. 2nd ed. Burlington, USA: Butterworth-Heinemann. 2007. ISBN: 978-0-7506-6927-6
14. 14. User’s Guide. Control System ToolboxTM. Revised for Version 9.7 (Release 2014a). The MathWorks, Inc; 2014
15. 15. Kisabo AB. Expendable Launch Vehicle Flight Control-Design & Simulation with Matlab/Simulink. Germany: Lambert Academic Publishing (LAP). ISBN: 978-3-8443-2729-8; 2011

Written By

Aliyu Bhar Kisabo and Aliyu Funmilayo Adebimpe

Reviewed: October 29th, 2018 Published: June 5th, 2019