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

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


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

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.

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: where F A xb , F A yb , F A zb ¼ components of aerodynamic force vector F A expressed in the body.
coordinate system, N.
F g xb , F g yb , F g zb ¼ components of gravitational force vector F g expressed in body coordinate system, N.
F p xb , F p yb , F p zb ¼ components of thrust vector F p 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 V x expressed in the body coordinate system, m/s. _ u, _ v, _ w ¼ components of linear acceleration expressed in body coordinate system, m/s 2 .
The aerodynamic forces have the following components: where C A = aerodynamic axial force coefficient, dimensionless.
C N = aerodynamic normal force coefficient, dimensionless.
C Ny = coefficient corresponding to component of normal force on yb-axis , dimensionless C Nz = coefficient corresponding to component of normal force on zb-axis , dimensionless S = aerodynamic reference area, m 2 .
V M = magnitude of velocity vector of the center of mass of the rocket, m/s. r = atmospheric density, kg/m 3 .
The propulsive forces in (1), as depicted in (Figure 1), are computed [8] as follows: where F p = magnitude of total instantaneous thrust force vector, N.
F pxb , F pyb , F pzb = components of thrust vector F p expressed in body coordinate system, N. γ 1 = angle measured from xb-axis to projection of thrust vector F p on x b y b -plane, rad (deg). γ 2 = angle measured from projection of thrust vector F p on x b y b -plane to the thrust vector F p , rad.
l p = 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, F p will be computed as given in (3). Here we assume that F p is acting along the line of symmetry of the rocket because the nozzle is fixed (fin control). Hence, the magnitude of the thrust force F p is calculated by where A e = rocket nozzle exit area, m 2 .
F pref = reference thrust force, N. P a = ambient atmospheric pressure, Pa.
The gravitational forces in (1) are computed as follows: where F gxe , F gye, F gze = components of gravitational force vector F g expressed in earth coordinate system, N. g = acceleration due to gravity, m/s 2 . m = instantaneous mass of rocket, kg.
The dependence of the acceleration due to gravity on the altitude of the rocket is given by where g = acceleration due to gravity, m/s 2 .
g 0 = acceleration due to gravity at earth surface (nominally 9.8 m/s 2 ), m/s 2 . h = altitude above sea level, m.
R e = radius of the earth, m.
The gravitational force expressed in body coordinates is computed by multiplying (5) by matrix (7): 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).
[T e/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 Hence, considering (5) we can write the following: The terms F gxb , F gyb , and F gzb are the components of the gravitational force substituted into (1).
The mass in (1) is given below: where F pref = reference thrust force, N.

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 L A , M A , N A = components of aerodynamic moment vector M A expressed in body coordinate system (roll, pitch, and yaw, respectively), Nm.
L p , M p , N p = components of propulsion moment vector M p expressed in body coordinate system (roll, pitch, and yaw, respectively), Nm. I x , I y , I z = components of inertia (diagonal elements of inertia matrix when products of inertia are zero), kg-m 2 .
where C l = aerodynamic roll moment coefficient about center of mass, dimensionless.
C m = aerodynamic pitch moment coefficient about center of mass, dimensionless.
C n = aerodynamic yaw moment coefficient about center of mass, dimensionless. d = aerodynamic reference length of body, m.
The aerodynamic moment coefficients are of the form where C lp = roll damping derivative relative to roll rate p, rad À1 (deg À1 ).
C lδ = slope of curve formed by roll moment coefficient C 1 versus control surface deflection, rad À1 (deg À1 ).
C mref = pitching moment coefficient about reference moment station, dimensionless.
C Ny = coefficient corresponding to component of normal force on yb-axis, dimensionless.
C Nz = coefficient corresponding to component of normal force on zb-axis, dimensionless.
C nref = yawing moment coefficient about reference moment station, dimensionless.
x cm = instantaneous distance from rocket nose to center of mass, 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 where C mref = 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. C m versus angle of attack, α rad À1 (deg À1 ) slope of tune formed.
C mδ = slope of tune formed by pitch moment coefficient C m versus control surface deflection for pitch, δ p rad À1 (deg À1 ).
C nβ = slope of curve formed by yawing moment coefficient C n versus angle of sideslip β, rad À1 (deg À1 ).
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 α 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.
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.

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

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 o la (state vectors) and o lo , for the lateral and longitudinal missile dynamics as 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 c la and c lo , described in (22). The control variables (input variables) describe the trimmed pilot control and the aircraft attitude: 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): In (22), J i is the Jacobian matrix evaluated near control input c i . Its entries are first-order partial derivatives and represent the effect of changes in each control input on acceleration.

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 firstorder differential equations and an algebraic output equation as [12]: x n ¼ f n x 1 ; …; x n ; u ð Þ y ¼ h x 1 ; …; x n ; u ð Þ The column vector x = [x 1 ,…x n ] 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 ð Þ ¼ Concisely, (26) is written as 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 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: 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 _