The Online Journal on Computer Science and Information Technology (OJCSIT) Vol. (1) No.

(3)

Design of Missile Control System

Using Model Predictive Control
M. Fawzy(*), M. A. S. Aboelela(+), O. Abd El Rhman(*), H. T. Dorrah(+)
(*)
Egyptian Armed Forces
(+)
Cairo University, Faculty of Engineering,
Electric Power and Machines Dept., Giza, Egypt.

Abstract-The goal of this paper is to control the II. MATHEMATICAL MODEL OF THE MISSILE
trajectory of the flight path of six degree of freedom flying
body model using Model predictive control (MPC) The model constitutes the six degree of freedom (6-DOF)
controller. MPC controller with constraints will be equations that break down into those describing kinematics,
developed and able to compensate for constraints that dynamics (aerodynamics, thrust, and gravity), command
represent physical limits of actuators in pitch and yaw guidance generation systems, and autopilot (electronics,
angles. The design of MPC controller with linear system instruments, and actuators). The input to this model is launch
for the six degree of freedom flying body is described. conditions, target motion, and target trajectory
MPC controllers are computationally intensive because an characterization, while the outputs are the missile flight data
on-line optimization problem is formed and solved at each (speed, acceleration, range, etc.) during engagement.
control cycle. The basic frames needed for subsequent analytical
developments are the ground, body and velocity coordinate
Keywords-six degree of freedom missile model, systems. The origins of these coordinate systems are the
Linearization model, Model predictive control. missile center of gravity (c.g). In the ground coordinate
system, the Xg-Zg lie in the horizontal plane and the Yg axis
I. INTRODUCTION completes a standard right-handed system and points up
vertically. In the body coordinate system, the positive Xb axis
In recent years, the requirements for the quality of coincides with the missile's center line and it is designated as
automatic control increased significantly due to increased
roll-axis. The positive Zb axis is to the right of the Xb axis in
complexity of plants and sharper specifications of product.
This paper will address the design of optimal variable the horizontal plane and it is designated as the pitch axis. The
structure controllers applied to a six degree of freedom positive Yb axis points upward and it is designed as the yaw
missile model. The six degree of freedom missile model is the axis. The body axis system is fixed with respect to the missile
solution to obtain a detailed accurate data about the missile and moves with the missile. In the velocity coordinate system,
trajectory. Linear model of the investigated systems will be Xv coincides with direction of missile velocity, (Vm) which
considered. The linearization will be obtained during two related to the directions of missile flight. The axis Zv
phases that are boost phase and sustain phase so the completes a standard right-handed system, [6, 11].
controllers are designed for two linear time invariant LTI The pitch plane is X-Y plane, the yaw plane is X-Z plane,
models. The paper objectives are and the roll plane is Y-Z plane. The ground coordinate system
To develop a general flexible sophisticated mathematical and body coordinate system are related to each other through
model of flight trajectory simulation for a hypothetical Eulers angles (, , ). The ground coordinate system and
anti tank missile, which can be used as a base line velocity coordinate system are related to each other through
algorithm contributing for design, analysis, and the (, ) angles. In addition, the velocity coordinate system
development of such a system and implement this model is related to the body frame through the angle of attack () in
using Simulink to facilitate the design of its control the pitch plane and angle of attack () in the yaw plane
system (sideslip angle). The angles between different coordinate
Developing control system, by using MPC control systems are shown in Fig. 1, [6, 11].
techniques The relation between the body and the velocity coordinate
systems can be given as follows
This paper is organized as follows. Section 2 reviews X b cos() cos( ) cos() sin ( ) sin () X v
mathematical model of six degree of freedom missile Y = sin ( ) cos( ) 0 Yv (1)
equations and linearization model is represented. Section 3 b
gives MPC controller design. Section 4 presents control Z b sin () cos( ) sin () sin ( ) cos() Z v
applications and results. Finally, conclusions are outlined in The body and velocity axes system as well as forces,
section 5. moments and other quantities are shown in Fig. 2.

Reference Number: W11-0028 64

The Online Journal on Computer Science and Information Technology (OJCSIT) Vol. (1) No. (3)

moments acting on missile in body coordinate system; x, y,

z are angular velocity in body coordinate system; Ix, Iy, Iz
are moments of inertia in body coordinate system; X is
missile range; Y is missile altitude; Z is horizontal
displacement of the missile; and m is missile mass. The forces
and the moments acting on missile are due to thrust, gravity
and aerodynamic forces and moments are given as following,
[4, 6, 11].
(
Fx = T cos( ) cos )
( ( ))
QS C x 0 + C x 2 + 2 mg sin ()
(14)

Fy = T sin ( ) + QS C y mg cos() (15)

Figure (1): The angles between different coordinate systems

( )
Fz = T cos( ) sin QS C z (16)
x D
M x = DQS m x 0 (17)
2Vm
( )
M y = T cos( ) sin X cg
y D (18)
+ DQS m y + m y0
Vm
D
M z = T sin ( )X cg + DQS m z + m z 0 z (19)
Vm
In these equations Cx, Cx0, Cy, Cz are aerodynamic force
coefficient; mx0, my, my0, mz, mz0 are aerodynamic
moment coefficients; D is the diameter of maximum cross
section area of body; S is the reference area; Q is the dynamic
Figure(2): Forces, moments and other quantities pressure; is the nozzle deflection angle in the pitch plane;
is the nozzle deflection angle in the yaw plane; T is the
There are 6 dynamic equations (3 for translational motion
and 3 for rotational motion) and 6 kinematic equations (3 for thrust force; Xcg is the distance between center of gravity (cg)
translational motion and 3for rotational motion) for a missile and the nozzle; and g is acceleration due to gravity and is
with six degrees of freedom. The equations are somewhat taken to be constant 9.81 m/sec2.
simpler, if the mass is constant. The missile 6DOF equations The linearized model takes the following form
in velocity coordinate system are given as following, [6] x= Ax + Bu (20)
Fx = m Vm (2) y = Cx + Du (21)
where
F = mV
y m (3)
[
x = Vm X Y Z x y z ]
Fz = m Vm cos()
u = [ ]
(4)
( )
x I y I z z y
M x = I x (5) t

y (I y I x ) z x
M y = I y (6) y = [ ]T
M z = Izz (I x I y ) x y (7) Where matrices A, B, C and D are matrix coefficient of LTI
system; x is state vector; u is input vector; y is output vector.
= V cos() cos( )
X (8) Missile solid propellant thrust will be divided into two
m
= V sin ()
Y (9)
phases, first phase is Boost phase that will take about 5.8 sec
m of total flight time (0 t < 5.8 sec) and thrust force T = Tmax.
Z = V cos()sin ( )

m (10) The second phase is Sustain phase that will start after boost
( )
= y cos( ) z sin ( ) cos() (11) region until the impact with target (5.8 t < 25 sec) and
thrust force T = Tmin. Therefore, we will discuss the
= y sin ( ) + z cos( ) (12) linearization of missile motion equation at boost and sustain
( )
= x tan () y cos( ) z sin ( ) = x sin () (13) phases.
A linear time-invariant (LTI) model is implemented in a
In these equations, Fx, Fy, Fz are component of forces acting boost phase around the operating point at t = 0. A linear time-
on missile in velocity coordinate system; Mx, My, Mz are invariant (LTI) model is implemented in sustain phase around

Reference Number: W11-0028 65

The Online Journal on Computer Science and Information Technology (OJCSIT) Vol. (1) No. (3)

the operating point at t = 5.8 sec. Open loop LTI model of known values up to instant k (past input and output) and
missile motion equation is represented in Fig.3, [2, 9]. on the future control signals u (k + i 1 k ) for i = 1 M
(M is control horizon for future moves and 1 M P)
III. MODEL PREDICTIVE CONTROL (MPC) DESIGN
which are those to be sent to the system and calculated.
The prediction future output can be explained briefly in
Model Predictive Control (MPC) is an advanced method of
[10, 12].
2. The set of future control signals u (k + i 1 k ) is calculated
process control that has been in use in the process industries
since the 1980s. MPC is a control strategy that is suitable for
optimizing the performance of constrained systems. by optimizing a determined criterion to keep the process as
Constrains are present in all control systems due to actuators close as possible to the reference trajectory r(k+i) that can
physical limits, boundaries of safe operation and lower limits be the set point itself or a close approximation of it. This
for product quality. The MPC uses the same powerful linear criterion usually takes the form of a quadratic function of
dynamic modeling that employ transfer functions, state-space the errors between the predicted output signal and the
matrices, or a combination of the two. MPC systems rely on predicted reference trajectory. The control effort is
the idea of generating values for process inputs as solutions of included in the objective function in most cases. An
an on-line (real-time) optimization problem. This problem is explicit solution can be obtained if the criterion is
based on a process model and process measurements, [5, 10, quadratic, the model is linear, and there are no constraints;
12]. otherwise an iterative optimization method has to be used.
Some assumptions about the structure of the future control
law are also made in some cases, such as that it will be
constant from a given instant.

Figure (3): Open loop LTI model of six degree-of-freedom Figure (4): Block diagram of SISO plant MPC controller
missile equation
3. The control signal u (k k ) is sent to the process whilst the
Fig. 4 shows block diagram of single input single output
(SISO) plant MPC controller. The main objective is to hold a next control signals calculated are rejected, because at the
single unmeasured output yu at a reference value (or set next sampling instant y(k+1) is already known and step 1
point), r, by adjusting a single manipulated variable (or is repeated with this new value and all the sequences are
actuator) u. The SISO plant actually has multiple inputs brought up to date. Thus the u (k + 1 k ) is calculated, which
in principle will be different from the u (k + 1 k ) because of
(manipulated variable input u, measured disturbance v and
unmeasured disturbance d). The controller receives the
measured disturbance v directly. This allows the controller to the new information available, using the receding horizon
compensate for measured disturbance impact on yu concept.
immediately rather than waiting until the effect appears in the
measured output ym. This is called feed forward control. MPC Fig. 5 shows the state of a hypothetical SISO MPC system
design always provides feedback compensation for that has been operating for many sampling instants. Integer k
unmeasured disturbances and feed forward compensation for represents the current sampling instant. The current measured
any measured disturbance, The MPC design removes the output, y(k), and previous measurements y(k1), y(k2), ...
estimated noise z component of the measurement (filtering). are known and are the filled circles in Fig. 5 (a), [7, 9]. Fig. 5
[9]. (b) shows the controllers previous moves, u(k4), u(k3), ,
The methodology of all the controllers belonging to the u(k-1) as filled circles. As is usually the case, a zero-order
MPC family is characterized by the following strategy, [5]: hold receives each move from the controller and holds it until
1. The future outputs for a determined prediction horizon P the next sampling instant, causing the step-wise variations
(future sampling periods P 1) are predicted at each shown in Fig. 5 (b).
sampling instant k using the process model. These For the basic formulation of predictive control, we shall
prediction outputs y(k + i k ) for i = 1 . P depend on the
assume that the plant model is linear, that the objective

Reference Number: W11-0028 66

The Online Journal on Computer Science and Information Technology (OJCSIT) Vol. (1) No. (3)

function (cost function) is quadratic, and that constraints are Where:

in the form of linear inequalities. Furthermore, we shall Q(i), R(i) are nonnegative optimization weight
assume that the cost function does not penalize particular coefficients.
values of the input vector u(k), but only changes of the input
vector, u (k ) = u (k ) u (k 1) . This formulation coincides 2
x Q = x T Qx
with that used in the majority of the predictive control
literature. To make the useful formulation, we shall not In some methods the second term, which considers the
assume the state variables can be measured, but that is control effort, is not taken into account, whilst in others the
obtained an estimate x (k k ) of the state x(k). Signals values of the control signal (not its increment) also appear
u (k + i 1 k ) will denote a future value of the input u. Signals
directly.
In practice all process are subject to constraints. The
x (k + i k ) and y(k + i k ) will denote the predictions, made at actuators have a limited field of action and a determined slew
rate, as is the case of the valves, limited by the positions of
time k, of the variables x and y at time k+i, [5].
totally open or closed and by the response rate. Constructive
reasons, safety or environmental ones, or even the sensor
scopes themselves can cause limits in the process variables
such as levels in tanks, flows in piping, or maximum
temperatures and pressures; moreover, the operational
conditions are normally defined by the intersection of certain
constraints for basically economic reasons, so that the control
system will operate close to the boundaries. All of this makes
the introduction of constraints in the function to be minimized
necessary. Normally, bounds in the amplitude and limits in
the output will be considered, [3, 12].
By adding constraints to the objective function, the
minimization becomes more complex, so that the solution
cannot be obtained explicitly as in the unconstrained case.
Constraints may be either hard or soft. A hard constraint must
not be violated. Unfortunately, under some conditions a
constraint violation might be unavoidable (e.g., an
unexpected, large disturbance), and a realistic controller must
allow for this. The MPC does so by softening each constraint,
making a violation mathematically acceptable, though
discouraged. The designer may specify the degree of softness
in each case, making selected constraints less likely to be
Figure (5): MPC controller at the kth sampling instant violated than others. Briefly, you specify a tolerance band for
each constraint. If the tolerance band is zero, the constraint is
A cost function J penalizes deviations of the predicted hard (no violation allowed). Increasing the tolerance band
controlled outputs y(k + i k ) from a reference trajectory softens the constraint. The tolerance band is not a limit on the
constraint violation, however. (If it were, you would still have
r(k+i). The reference trajectory may depend on measurements
a hard constraint.) You need to view it relative to other
made up to time k, in particular, its initial point is the output
constraints, [9].
measurement y(k). The reference trajectory may also be a
fixed set point, or some other predetermined trajectory. The
VI. CONTROL APPLICATION AND RESULTS
general expression for such an objective function (cost
function) will be, [5, 10, 12]:
In this section, the autonomous flight of six degree of
J = y(k + i k ) r (k + i )
P 2 freedom flying body is simulated. The goal is to control the
Q (i ) trajectory of the flight path of six degree of freedom flying
i =1
(22) body model using MPC controller. The design of MPC
+ u (k + i 1 k ) r (k + i )
M 2 controller with LTI system for the six degree of freedom
R (i ) flying body is described. This design has been implemented
i =1
u min u (k + i 1 k ) u max
in a simulation environment under Matlabs toolbox
i = 1, , M
u (k + i 1 k ) u max
Simulink.
u i = 1, , M A linear time-invariant (LTI) model is implemented around
subject to min (23)
y min y(k + i k ) y max i = 1, , P the operating point at =. The state and input at this point are
u (k + i 1 k ) = 0 i = M + 1, , P described as following

Reference Number: W11-0028 67

The Online Journal on Computer Science and Information Technology (OJCSIT) Vol. (1) No. (3)

x 0 = [8.6 0.6108 0 0 0 0 0 0.6108 0 0 0 0] problem of controller tuning can be solved by an application

of the particle swarm optimization (PSO) algorithm for
u 0 = [0 0 0] optimization on a three-dimensional solution space, each
particle having a three-dimensional position and velocity
The nozzle deflection angle in pitch plane ( ) and yaw vector. The initial positions of the ith particles of the swarm
can be represented by a three dimensional vector, and then the
plane ( ) are limited with 28.5 (0.5rad).
initial values are randomly generated based on the extreme
The pitch demand programmer is an exponential command values. PSO algorithm and its parameters chosen can be
and is described as explained briefly in [8, 13].
t p
p = p0 A Be
The position vector of the best particle gives the optimized
(24)
parameter of integer PID controller as following, [8]
where p0 is the missile-launching angle with respect to the The PID controller parameters for pitch angle during
horizon; A, B are vertical position angles depending on boost phase are kp = 20.432, ki = 4.1353, kd = 1.241.
The PID controller parameters for pitch angle during
target position. For our simulation p0 = 35 ;
sustain phase are kp = 30.2494, ki = 20.6635, kd = 5.4857.
A = B = 30; p = 2.1788 sec . . The PID controller parameters for yaw angle during boost
phase are kp = 20.432, ki = 4.1353, kd = 1.241.
The yaw demand programmer is an exponential command
The PID controller parameters for yaw angle during
and is described as
sustain phase are kp = 30.2494, ki = 20.6635, kd =
t
p = s 1 e (25) 5.4857.

where s is horizontal position angle depending on target The pitch error is the difference between pitch demand
program (pitch reference trajectory) and pitch angle response.
position. For our simulation s = 5; = 0.2 sec . Fig.8 represents the pitch error comparison with PID and
MPC. The pitch error with MPC has high overshoot and high
To design an MPC controller the MPC toolbox in Matlabs oscillation at starts of boost and sustain phases (at t = 0 and
toolbox Simulink has been used. The controller design t = 5.8sec) and also it has small steady state error during
requires a LTI model of the plant that is to be controlled. sustain phase However, for pitch error with PID controller
Multi input multi output MPC controller is designed where has small overshoot and there is no oscillation.
one MPC controller is used during boost and sustain phases. The yaw error is the difference between yaw demand
The rate at which MPC operates is 1/NTS, where TS is program (yaw reference trajectory) and yaw angle response.
control interval (sampling period), N is the number of Fig.8 represents the yaw error with PID and MPC. The yaw
controls in the control history that are applied to the plant. error of nonlinear system with MPC has high overshoot at
N = 1 is chosen since this is the value N typically takes. The start of boost (at t = 0) and it has high oscillation at start of
value of TS is important since it is the length of each sustain phase (t = 5.8sec), also the steady state error increases
prediction step and the duration each control input is held during sustain phase. However, for yaw error with PID
constant. The method for choosing TS for this problem is controller has small overshoot.
based on tracking performance. After further tests TS = 0.01 is
chosen since this value of TS gives the best tracking V. CONCLUSION
performance to a sequence of pitch and yaw demand
generator programs. [1] The design of PID controllers gave the best response for
Selecting the prediction horizon P was also affected by the pitch and yaw angles where there is no oscillation (chattering)
controller. To keep the controller simple, the prediction and has small overshoot. The parameters optimization of PID
horizon P and control horizon M were set equal to each other. controllers based on PSO method was highly effective.
After further tests a value of P = M = 1 was chosen since According to optimization target, the PSO method could
these value of P, M give the best tracking performance to a search the best global solution for PID controllers parameters
sequence of pitch and yaw demand generator program. [1] and guarantee the objective solution space in defined search
Optimization parameters (Q, R) start with identity matrices; space.
the values were changed through trial and error to improve The design of MPC gave response less quality than that
the tracking performance to a sequence of pitch and yaw was given from PID controller but acceptable responses.
demand generator program [1]. After further tests a values of However, MPC controller can be used to control a great
Q, R was chosen as following variety of processes (one MPC controller was used instead of
four PID controllers). MPC controller is a limited knowledge
3.1899 0
Q= R=0 of control because its tuning is relatively easy. MPC
3.1899
, (26)
0 controller is simple to the treatment of constraints and relies
The PID controller has three unknown parameter kp, ki and on the idea of generating values for process inputs as
kd that are required to be designed. Hence, the present solutions of an on-line (real-time) optimization problem.

Reference Number: W11-0028 68

The Online Journal on Computer Science and Information Technology (OJCSIT) Vol. (1) No. (3)

Figure (6): Pitch and yaw angles with MPC controller versus time

Figure (7): Pitch and yaw angles with PID controller versus time

Figure (8): Pitch error and yaw error comparisons with PID and MPC

Reference Number: W11-0028 69

The Online Journal on Computer Science and Information Technology (OJCSIT) Vol. (1) No. (3)

REFERENCES [8] D. Maiti, Acharya, A., Chakraborty, M. and Konar, A.,

[1] A. Alaniz, Model Predictive Control with Application Tuning PID and PID Controllers using the Integral
to Real Time Hardware and a Guided Parafoil, Master Time Absolute Error Criterion, IEEE, 2008.
degree of Science in Aeronautics and Astronautics, [9] MATLAB 9.0 - User's Guide, the Math Works Inc.,
Massachusetts Institute of Technology, June 2004 Natick, MA, USA, 2010.
[2] T. Ashish, Modern Control Design with MATLAB and [10] P. E. Orukpe, Basics of Model Predictive Control,
SIMULINK Indian Institute of Technology, Kanpur, Master degree of Science in Control Engineering,
India, John Wiley & Sons, 2002. Imperial College, London, April 2005.
[3] A. Bemporad, Morari, M., Dua, V. and Pistikopoulos, E. [11] G. M. Siouris, Missile Guidance and Control Systems,
N., The Explicit Solution of Model Predictive Control Springer-Verlag, New York, USA, 1st edition, March
via Multi parametric Quadratic Programming, 2004.
American Control Conference, Chicago, Illinois, June [12] A. A. Tyagunov, High-Performance Model Predictive
2000. Control for Process Industry, Ph.D. degree of Science
[4] J. H. Blakelock, Automatic Control of Aircraft and in Control system, Faculty of Electrical Engineering,
Missiles, John Wiley & Sons, Inc., USA, 2nd edition, Eindhoven University of Technology, June 2004.
1991. [13] J. Yi Cao and Gang Cao, B., Design of Fractional
[5] E. F. Camacho and Bordons, C., Model Predictive Order Controller Based on Particle Swarm
Control, Springer-Verlag, London, 2nd edition, 2004. Optimization, International Journal of Control,
[6] P. Garnell, Guided weapon control systems, Pergamon Automation and Systems, Vol. 4, No. 6, pp. 775-781,
Press, Oxford, New York, 2nd edition, 1980. December 2006.
[7] M. Lazar, Model Predictive Control of Hybrid Systems:
Stability and Robustness, Ph.D. degree of Science in
Control system, Faculty of Electrical Engineering,
Eindhoven University of Technology, September 2006.

Reference Number: W11-0028 70