Open access peer-reviewed chapter

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

By Aliyu Bhar Kisabo and Aliyu Funmilayo Adebimpe

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

DOI: 10.5772/intechopen.82292

## 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 .

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 . Especially, in the field of control system design, a mathematical model in transfer function might not be ideal for optimal controldesign. 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.”

## 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  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-forceor 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̇=FAxb+FPxb+Fgxbmqwrv,m/s2v̇=FAyb+FPyb+Fgybmrupw,m/s2ẇ=FAzb+FPzb+Fgzbmpvqu,m/s2E1

where

FAxb,FAyb,FAzb=components of aerodynamic force vector FAexpressed in the body.

coordinate system, N.

Fgxb,Fgyb,Fgzb=components of gravitational force vector Fgexpressed in body coordinate system, N.

Fpxb,Fpyb,Fpzb=components of thrust vector Fpexpressed 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:

FAxb=0.5ρVM2CASFAyb=0.5ρVM2CNySFAzb=0.5ρVM2CNzSE2

where

CA = aerodynamic axial force coefficient, dimensionless.

CN = aerodynamic normal force coefficient, dimensionless.

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

CNy=CNvv2+w2,dimensionless

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

CNz=CN1v2+w2,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  as follows:

Fpxb=Fpcosγ2cosγ1,NFpyb=Fpcosγ2sinγ1,NFpzb=Fpsinγ2,NE3

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 Fpon xb yb-plane, rad (deg).

γ2= angle measured from projection of thrust vector Fpon 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 Fpis calculated by

Fp=Fpref+prefpaAe,NE4

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:

Fgxe=0,NFgye=0,NFgze=mg,NE5

where

Fgxe, Fgye, Fgze = components of gravitational force vector Fgexpressed 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=g0Re2Re+h2,m/s2E6

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):

Te/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θ,dimensionlessE7

where

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

s = sinefunction (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 vexpressed in the body coordinate system can be transformed to the earth coordinate system by the matrix equation

ve=Te/bvbE8

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

FgxbFgybFgzb=Tb/eFgxeFgyeFgze,NE9

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

The massin (1) is given below:

m=m01Isp0tFprefdt,kgE10

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

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).

LA=0.5ρVM2ClSd,NmMA=0.5ρVM2CmSd,NmNA=0.5ρVM2CnSd,NmE12

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

Cl=Clδδr+d2VMClppCm=CmrefCNzxcmxrefd+d2VMCmq+CmαqCn=Cnref+CNyxcmxrefd+d2VMCnr+CnβrdimensionlessE13

where

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

C = slope of curve formed by roll moment coefficient C1versus 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).

Cmα̇= 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.

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

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

Cnβ̇= 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

Cmref=Cmαα+CmδδηCnref=Cnββ+Cnδδζ,dimensionlessE14

where

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

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

Cmδ=slope of tune formed by pitch moment coefficient Cmversus control surface deflection for pitch, δprad−1 (deg−1).

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

Cnδ=slope of curve formed by yaw moment coefficient Cnversus effective control surface deflection for yaw, δyrad−1 (deg−1).

α= angle of attack, rad (deg).

β= angle of sideslip, rad (deg).

δη = δ1= δ3effective 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

where ∅ is aerodynamic roll angle, rad (deg), and αtis 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®  and Flexible Structures Simulator (FSS) . 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/NVariableDescription
1CNNormal force coefficient (body axis)
2CLLift coefficient (wind axis)
3CMPitching moment coefficient (body axis)
4XcpCenter of pressure in calibers from moment reference center
5CAAxial force coefficient (body axis)
6CDDrag coefficient (wind axis)
7CYSide force coefficient (body axis)
8CnYawing moment coefficient (body axis)
9ClRolling moment coefficient (body axis)
10CNormal force coefficient derivative with angle of attack
11CPitching moment coefficient derivative with angle of attack
12CSide force coefficient derivative with sideslip angle
13CYawing moment coefficient derivative with sideslip angle (body axis)
14CRolling moment coefficient derivative with sideslip angle (body axis)
15CMqPitching moment coefficient derivative with pitch rate
16CNqNormal force coefficient derivative with pitch rate
17CAqAxial force coefficient derivative with pitch rate
18CMα̇Pitching moment derivative with rate of change of angle of attack
19CNα̇Normal force derivative with rate of change of angle of attack
20ClpRolling moment coefficient derivative with roll rate
21CnpYawing moment coefficient derivative with roll rate
22CYpSide force coefficient derivative with roll rate
23ClrRolling moment coefficient derivative with yaw rate
24CnrYawing moment coefficient derivative with yaw rate
25CYrSide 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.,

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̇=FAxb+FPxb+Fgxbmqwrv,m/s2v̇=FAyb+FPyb+Fgybmrupw,m/s2ẇ=FAzb+FPzb+Fgzbmpvqu,m/s2

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

ola=v̇ṗṙβT.E18
olo=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 claand clo, described in (22). The control variables (input variables) describe the trimmed pilot control and the aircraft attitude:

cla=δaδrϕψT,E20
clo=δeτθT,E21

Finally, the 12 states of the 6DoF equation of motion must be initialized with the initial state conditions. In MATLAB®, the trimcommand 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):

oi+1oi=Jici+1ci.E22

In (22), Jiis 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:

Ji=occila=u̇δau̇δru̇ϕu̇ψv̇δav̇δrv̇ϕv̇ψẇδaẇδrẇϕẇψṗδaṗδrṗϕṗψq̇δaq̇δrq̇ϕq̇ψṙδaṙδrṙϕṙψβδaβδrβϕβψci.E23
Ji=occilo=u̇δeu̇τu̇θv̇δev̇τv̇θẇδeẇτẇθṗδeṗτṗθq̇δeq̇τq̇θṙδeṙτṙθβδeβτβθci.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 nfirst-order differential equations and an algebraic output equation as :

ẋ1=f1x1xnuẋ2=f2x1xnuẋn=fnx1xnuy=hx1xnuE25

The column vector x = [x1,…xn]T is called the stateof the system. The scalars uand yare called the control input and the system output, respectively, denoting

fxu=f1x1xnuf2x1xnufnx1xnu.E26

Concisely, (26) is written as

ẋ=fxu,y=hxu.E27

where fand hare nonlinear functions of xand u; then, we say that the system is nonlinear. To linearize (26), we desire it to become

where Ais n x n, Bis n x 1, Cis 1 x n, and Dis 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 = [x1xn]Tobtained when u = u, noting that x=x=x,we define a coordinate transformation as follows:

Δx=Δx1Δxn=x1x1xnxn.

Further, denotingu=uu, y=yhxu. The new coordinates x,u,and yrepresent the variations of x, u, and yfrom 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  state is given by

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

where

A=fux,u=f1x1x1xnuf1xnx1xnufnx1x1xnufnxnx1xnu,
B=fux,u=f1ux1xnufnux1xnu,
C=hxx,u=hx1x1xnuhxnx1xnu,D=hux,u.

Note that the matrices A, B, C, and Dhave constant coefficients in that all partial derivatives are evaluated at the numerical values x1xnu.

### 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)):

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

v̇ṗṙϕ̇ψ̇=yvypyryϕyψlvlplrlϕlψnvnpnrnϕnψ0100000100vprϕψ+yξyζlξlζnξnζ0000ξζ.E31

All the variables in Amatrix of (31) are the lateral dimensionless aerodynamic stability derivatives with respect to the system state vectors. The variables in Bmatrix 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 (Fxand Fz), two of the momentum equations (Mx and Mz), and two of the Euler angles from the 6DoF equations as given in (17):

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

u̇ẇq̇θ̇=xuxwxqxθzuzwzqzθmumwmqmθ0010uwqθ+xηxτzηzτmηmτ00ητE33

In MATLAB® the linmod 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 xzplane in which the inertia cross-coupling in (xyand xzplanes) 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  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̇=01014.780500.01958100.85810.1256θqw+03.485820.42δη+014.780594.8557αwy=100θqwE34

where

A=01014.780500.01958100.85810.1256,B=03.485820.42,G=014.780594.8557C=100,D=0.

## 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.

## 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.

chapter PDF
Citations in RIS format
Citations in bibtex format

## More

© 2019 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution 3.0 License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

## How to cite and reference

### Cite this chapter Copy to clipboard

Aliyu Bhar Kisabo and Aliyu Funmilayo Adebimpe (June 5th 2019). State-Space Modeling of a Rocket for Optimal Control System Design, Ballistics, Charles Osheku, IntechOpen, DOI: 10.5772/intechopen.82292. Available from:

### Related Content

#### Ballistics

Edited by Charles Osheku

Next chapter

#### Discrete Element Modeling of a Projectile Impacting and Penetrating into Granular Systems

By Waseem Ghazi Alshanti

First chapter

#### Introductory Chapter: Laminations - Theory and Applications

By Charles Attah Osheku

We are IntechOpen, the world's leading publisher of Open Access books. Built by scientists, for scientists. Our readership spans scientists, professors, researchers, librarians, and students, as well as business professionals. We share our knowledge and peer-reveiwed research papers with libraries, scientific and engineering societies, and also work with corporate R&D departments and government entities.