# Autopilot Design Method for the Blended Missile Based on Model Predictive Control

Authors: Baoqing Yang; Yuyu Zhao
Journal: International Journal of Aerospace Engineering
 (2015)

---

## Abstract

This paper develops a novel autopilot design method for blended missiles with aerodynamic control surfaces and lateral jets. Firstly, the nonlinear model of blended missiles is reduced into a piecewise affine (PWA) model according to the aerodynamics properties. Secondly, based on the equivalence between the PWA model and mixed logical dynamical (MLD) model, the MLD model of blended missiles is proposed taking into account the on-off constraints of lateral pulse jets. Thirdly, a hybrid model predictive control (MPC) method is employed to design autopilot. Finally, simulation results under different conditions are presented to show the effectiveness of the proposed method, which demonstrate that control allocation between aerodynamic control surfaces and lateral jets is realized by adjusting the weighting matrix in an index function.


---

## Body

## 1. Introduction

The emergence of highly maneuverable targets has brought new challenges to guidance technology. The improvement of homing guidance performance against highly maneuverable targets in the future guided missiles requires the control system to have faster response and wider operation range [1]. In the meantime, as an important property of advanced missiles, multiple actuators are often employed to enhance maneuverability as well as interception probability. For example, both aerodynamic surfaces and reaction jets are employed in the control system of PAC-3.With a higher angle of attack, the missile dynamic model is highly nonlinear and the coupling effects [2] as well as the uncertainties in both aerodynamic parameters and reaction jet thrust are oblivious. Munson and Garbrick introduced amplification factors to describe the lateral jet interference effect [3]. Graham and Weinacht studied the interaction between the side jet and the external flow by a numerical method [4]. It is obvious that the autopilot designed using the linearized model around an operation point is usually unable to achieve satisfactory performance over a full flight envelope. Recently, many nonlinear control methods are proposed for the design of conventional aerofin autopilots. For instance, the gain-scheduled approach was proposed in [5], where linear parameter varying transformations were adopted.On the other hand, the autopilot design for blended missile has not been much reported, which is more complicated than that of conventional aerofin autopilot due to the heterogeneous actuation [1]. Hirokawa et al. [1] designed an autopilot for the case with aerodynamic surfaces and reaction jet using the coefficient diagram method (CDM). The feasibility of the autopilot based on variable structure was discussed in [6], where a simple blending strategy is investigated: aerodynamic control is used at low angles of attack while reaction jet control is used when the missile moves beyond the stall region. Optimal control and static control allocation were combined in [7] to address the issue of dual control missile while a dynamic control allocation method was presented in [8].In addition, for the autopilot design of blended missiles (such as PAC-3), two aspects of control input constraints should be taken into account: the saturation constraint on aerodynamic surfaces and the finite set constraint on reaction jets. However, they are not simultaneously considered in the aforementioned work. Besides, the hybrid properties of control inputs (continuous aerodynamic surfaces and on-off reaction jets) are often neglected as in [6–8]. Motivated by these facts, this paper attempts to design an autopilot using explicit hybrid MPC for blended missiles, by noting that MPC is a promising methodology for the control problem of constrained uncertain systems [9] and that the computational burden of on-line optimization is effectively reduced by using explicit MPC instead of traditional MPC [10].The remainder of this paper is as follows. Section2 gives a mathematical model of blended missile including the configuration of reaction jets. In Section 3, the piecewise affine model of blended missile is established, followed by an MLD model which is obtained based on the equivalence of piecewise affine model and mixed logical dynamical model. In Section 4, a hybrid MPC based method for autopilot design is proposed and an explicit control law is constructed. In Section 5, the effectiveness of the proposed method is verified by simulation cases under different conditions. Finally, several concluding remarks are given in Section 6.


## 2. Mathematical Model of the Missile

The plane reference coordinate systemO
U
V
W, the body coordinate system O
x
b
y
b
z
b, the trajectory coordinate system O
x
w
y
w
z
w, and velocity coordinate system O
x
v
y
v
z
v are involved in this paper. Figure 1 shows a missile with some key variables and identified axes. The axes z
w, z
v are not given, whose directions can be determined by the right hand rule.Figure 1
Key coordinate systems.


### 2.1. Missile Dynamic Model

The nonlinear motion equations are given by(1)
m
V
˙
=
P
cos
⁡
α
cos
⁡
β
-
X
a
-
m
g
sin
θ
+
F
x
w
a
,
m
V
θ
˙
=
P
sin
α
cos
⁡
γ
v
+
cos
⁡
α
sin
β
sin
γ
v
+
Y
a
cos
⁡
γ
v
-
m
g
cos
⁡
θ
+
F
y
w
a
,
m
V
cos
⁡
θ
ψ
˙
v
=
-
P
(
sin
α
sin
γ
v
-
cos
⁡
α
sin
β
cos
⁡
γ
v
)
-
Y
a
sin
γ
v
-
Z
a
cos
⁡
γ
v
-
F
z
w
a
,
(2)
J
x
p
˙
=
J
y
-
J
z
r
q
+
M
x
,
J
y
q
˙
=
J
z
-
J
x
p
r
+
M
y
,
J
z
r
˙
=
J
x
-
J
y
p
q
+
M
z
,
where F
x
w
a, F
y
w
a, F
z
w
a are lateral forces in the x, y, and z directions of trajectory coordinate system, respectively. For simplicity, we suppose that the body is symmetric about the x-axis; that is, J
x
y
=
J
y
z
=
J
z
x
=
0. In (2), each of the moments M
x, M
y, and M
z contains two components that are generated by aerodynamic surface and the lateral pulse jets, respectively:
(3)
M
x
=
M
e
x
+
M
x
b
a
,
M
y
=
M
e
y
+
M
y
b
a
,
M
z
=
M
e
z
+
M
z
b
a
,
where M
e
x, M
e
y, and M
e
z denote aerodynamic moment components, and M
x
b
a, M
y
b
a, and M
z
b
a are pulse jet moment components.


### 2.2. Lateral Jet Forces and Moments Model

As shown in Figure1, the lateral jet force is generated by 180 pulse jets located in front of the center of mass of PAC-3. These jets are divided into 10 rings and arranged in staggered positions (18 pulse jets are included in each ring). In each ring, these jets are uniformly distributed and the central angle between two neighboring jets is 20 degree. Use i (i
=
1,2
,
…
,
10) and j (j
=
1,2
,
…
,
18) to denote the ring’s label and the jet’s label in each ring, respectively. The distance between ring i and the center of mass is denoted by l
i, while spacing of adjacent rings is Δ
l. The layout scheme of pulse jets is shown in Figures 2 and 3. Assume the force generated by each individual jet is a constant F
m. In body coordinate system, the lateral force generated by the (
i
,
j
) pulse jet is given by
(4)
F
x
b
i
,
j
F
y
b
i
,
j
F
z
b
i
,
j
=
0
F
m
cos
⁡
2
j
-
i
*
18
π
-
F
m
sin
2
j
-
i
*
18
π
.
The corresponding moment is given by
(5)
M
x
b
i
,
j
M
y
b
i
,
j
M
z
b
i
,
j
=
0
F
m
l
i
sin
2
j
-
i
*
18
π
F
m
l
i
cos
⁡
2
j
-
i
*
18
π
,
where
(6)
i
*
=
2
,
i
is
odd
,
i
*
=
1
,
i
is
even
.The layout scheme of lateral pulse jets.
(a)
Odd-numbered rings
(b)
Even-numbered ringsFigure 3
The ring frames expansion of lateral pulse jets.For the situation where all pulse jets are fired at the same time, the total force and moment are given by(7)
F
x
b
F
y
b
F
z
b
=
0
∑
i
=
1
i
=
10
‍
∑
j
=
1
j
=
18
‍
F
y
b
i
,
j
∑
i
=
1
i
=
10
‍
∑
j
=
1
j
=
18
‍
F
z
b
i
,
j
M
x
b
M
y
b
M
z
b
=
0
-
∑
i
=
1
i
=
10
‍
∑
j
=
1
j
=
18
‍
F
z
b
i
,
j
l
i
∑
i
=
1
i
=
10
‍
∑
j
=
1
j
=
18
‍
F
y
b
i
,
j
l
i
.In order to avoid the coupling between the pitch moment and yaw moment, each ring is divided into four control regions: positive pitch, negative pitch, positive yaw, and negative yaw control region, as shown in Figure4.Schematic of control regions.
(a)
Odd-numbered rings
(b)
Even-numbered ringsThe autopilot design of PAC-3 is more complicated than that of other conventional missiles which are controlled only by aerodynamic surfaces, due to the hybrid property of control inputs and the on-off property of pulse jet (the pulse jet can be fired only one time). To deal with this problem, the work [7] proposed two-step design procedures: (1) in the first step, neglect the hybrid property (or on-off property) and design the expected force and moment signals, (2) in the other step, the fire logic is derived by solving the retaliation problem of these signals, with taking into account the hybrid property (or on-off property). Different from the above traditional procedures, a novel procedure will be presented in what follows, where only one step is included. In practical applications, only a small number of jets (in a certain of rings) are activated over a finite time interval. To make the idea of the following development clear, we here consider a simple but representative situation, where no more than two rings are allowed to be fired simultaneously and no more than two jets are activated in each fired ring. Meanwhile, it should be ensured that only odd rings or even rings are fired, and the jets are fired symmetrically about the corresponding symmetry axis of each control region.Take the positive pitch control region as an example. The forces provided by the jets(
i
,
1
), (
i
,
2
), (
i
,
3
), (
i
,
17
), and (
i
,
18
) in an odd ring are given by
(8)
F
o
=
F
m
F
m
cos
⁡
π
9
F
m
cos
⁡
2
π
9
F
m
cos
⁡
2
π
9
F
m
cos
⁡
π
9
T
.Similarly, the forces provided by the jets(
i
,
1
), (
i
,
2
), (
i
,
17
), and (
i
,
18
) in an even ring are given by
(9)
F
e
=
F
m
cos
⁡
π
18
F
m
cos
⁡
π
6
F
m
cos
⁡
π
6
F
m
cos
⁡
π
18
T
,
where F
o and F
e denote the forces associated with an odd ring and an even ring, respectively.UseF
i to denote the lateral force generated by ring i. Clearly, F
i should satisfy the condition
(10)
F
1
=
F
9
,
F
3
=
F
7
,
F
2
=
F
10
,
F
4
=
F
8
.When odd ring is fired, the lateral force and moment are given by(11)
F
y
b
∈
F
m
,
2
F
m
cos
⁡
π
9
,
2
F
m
cos
⁡
2
π
9
,
2
F
m
,
4
F
m
cos
⁡
π
9
,
4
F
m
cos
⁡
2
π
9
,
M
z
b
∈
F
m
l
5
,
2
F
m
l
5
cos
⁡
π
9
,
2
F
m
l
5
cos
⁡
2
π
9
,
2
F
m
l
5
,
4
F
m
l
5
cos
⁡
π
9
,
4
F
m
l
5
cos
⁡
2
π
9
.While even ring is fired, the lateral force and moment are given by(12)
F
y
b
∈
2
F
m
cos
⁡
π
18
,
2
F
m
cos
⁡
π
6
,
4
F
m
cos
⁡
π
18
,
4
F
m
cos
⁡
π
6
M
z
b
∈
2
F
m
l
6
cos
⁡
π
18
,
2
F
m
l
6
cos
⁡
π
6
,
4
F
m
l
6
cos
⁡
π
18
,
4
F
m
l
6
cos
⁡
π
6
.NotingΔ
l is quite small, we consider l
5
≃
l
6
=
l. In order to ensure that the jets fire efficiency, the jets resulting in small moment components along z
b axis are not activated. As a result, the sets of forces and moments are, respectively, given by
(13)
U
F
y
+
=
F
m
,
2
F
m
cos
⁡
π
9
,
2
F
m
cos
⁡
π
18
,
2
F
m
,
4
F
m
cos
⁡
π
9
,
4
F
m
cos
⁡
π
18
,
(14)
U
M
y
+
=
F
m
l
,
2
F
m
l
cos
⁡
π
9
,
2
F
m
l
cos
⁡
π
18
,
2
F
m
l
,
4
F
m
l
cos
⁡
π
9
,
4
F
m
l
cos
⁡
π
18
,
where F
y
b
∈
U
F
y
+ and M
z
b
∈
U
M
y
+. During each control period, a control moment belonging to U
M
y
+ will be used as the input.The mutual interference between high-speed jet stream and air leads to lateral jet interference effect. In order to take into account this interference, force and moment amplification factorsK
F
y, K
F
z, K
M
y, and K
M
z are introduced as in [3, 4]. Then, the resulting lateral forces and moments are
(15)
F
x
b
a
F
y
b
a
F
z
b
a
=
0
F
y
b
+
K
F
y
F
y
b
F
z
b
+
K
F
z
F
z
b
,
M
x
b
a
M
y
b
a
M
z
b
a
=
0
M
y
b
+
K
M
y
M
y
b
M
z
b
+
K
M
z
M
z
b
.Remark 1.
In fact, each pulse jet can be fired only once, so the location of the fired jet cannot provide force anymore. Based on this precondition, elements of setsU
F
y
+ and U
M
y
+ will be less and less over time. In this paper, quantity change of sets’ elements is not considered to simplify the problem.


### 2.3. Attitude Control Model

Some transformation and simplification are applied to the missile model for control design. It is assumed that the missile’s mass is of a constant value. Note that the goal is to establish the angle of attack and sideslip angle. The attitude control model is given by(16)
α
˙
=
r
+
q
sin
α
tan
β
-
Q
S
C
y
α
α
+
C
y
δ
z
δ
z
cos
⁡
α
m
V
cos
⁡
β
-
F
y
b
+
K
F
y
F
y
b
cos
⁡
α
m
V
cos
⁡
β
-
G
y
cos
⁡
α
m
V
cos
⁡
β
,
β
˙
=
q
cos
⁡
α
+
(
Q
S
(
C
z
β
β
+
C
z
δ
y
δ
y
)
+
F
z
b
+
K
F
z
F
z
b
)
cos
⁡
β
m
V
+
(
Q
S
(
C
y
α
α
+
C
y
δ
z
δ
z
)
+
F
y
b
+
K
F
y
F
y
b
)
sin
α
sin
β
m
V
+
G
z
cos
⁡
β
m
V
+
G
y
sin
α
sin
β
m
V
,
q
˙
=
M
y
b
J
y
+
K
M
y
M
y
b
J
y
+
Q
S
L
m
y
β
β
J
y
+
Q
S
L
m
y
δ
y
δ
y
J
y
+
Q
S
L
m
y
q
q
J
y
,
r
˙
=
M
z
b
J
z
+
K
M
z
M
z
b
J
z
+
Q
S
L
m
z
α
α
J
z
+
Q
S
L
m
z
δ
z
δ
z
J
z
+
Q
S
L
m
z
r
r
J
z
,
where m
y
β, m
y
δ
z, m
y
q, m
z
α, m
z
δ
z, and m
z
r are aerodynamic parameters.


## 2.1. Missile Dynamic Model

The nonlinear motion equations are given by(1)
m
V
˙
=
P
cos
⁡
α
cos
⁡
β
-
X
a
-
m
g
sin
θ
+
F
x
w
a
,
m
V
θ
˙
=
P
sin
α
cos
⁡
γ
v
+
cos
⁡
α
sin
β
sin
γ
v
+
Y
a
cos
⁡
γ
v
-
m
g
cos
⁡
θ
+
F
y
w
a
,
m
V
cos
⁡
θ
ψ
˙
v
=
-
P
(
sin
α
sin
γ
v
-
cos
⁡
α
sin
β
cos
⁡
γ
v
)
-
Y
a
sin
γ
v
-
Z
a
cos
⁡
γ
v
-
F
z
w
a
,
(2)
J
x
p
˙
=
J
y
-
J
z
r
q
+
M
x
,
J
y
q
˙
=
J
z
-
J
x
p
r
+
M
y
,
J
z
r
˙
=
J
x
-
J
y
p
q
+
M
z
,
where F
x
w
a, F
y
w
a, F
z
w
a are lateral forces in the x, y, and z directions of trajectory coordinate system, respectively. For simplicity, we suppose that the body is symmetric about the x-axis; that is, J
x
y
=
J
y
z
=
J
z
x
=
0. In (2), each of the moments M
x, M
y, and M
z contains two components that are generated by aerodynamic surface and the lateral pulse jets, respectively:
(3)
M
x
=
M
e
x
+
M
x
b
a
,
M
y
=
M
e
y
+
M
y
b
a
,
M
z
=
M
e
z
+
M
z
b
a
,
where M
e
x, M
e
y, and M
e
z denote aerodynamic moment components, and M
x
b
a, M
y
b
a, and M
z
b
a are pulse jet moment components.


## 2.2. Lateral Jet Forces and Moments Model

As shown in Figure1, the lateral jet force is generated by 180 pulse jets located in front of the center of mass of PAC-3. These jets are divided into 10 rings and arranged in staggered positions (18 pulse jets are included in each ring). In each ring, these jets are uniformly distributed and the central angle between two neighboring jets is 20 degree. Use i (i
=
1,2
,
…
,
10) and j (j
=
1,2
,
…
,
18) to denote the ring’s label and the jet’s label in each ring, respectively. The distance between ring i and the center of mass is denoted by l
i, while spacing of adjacent rings is Δ
l. The layout scheme of pulse jets is shown in Figures 2 and 3. Assume the force generated by each individual jet is a constant F
m. In body coordinate system, the lateral force generated by the (
i
,
j
) pulse jet is given by
(4)
F
x
b
i
,
j
F
y
b
i
,
j
F
z
b
i
,
j
=
0
F
m
cos
⁡
2
j
-
i
*
18
π
-
F
m
sin
2
j
-
i
*
18
π
.
The corresponding moment is given by
(5)
M
x
b
i
,
j
M
y
b
i
,
j
M
z
b
i
,
j
=
0
F
m
l
i
sin
2
j
-
i
*
18
π
F
m
l
i
cos
⁡
2
j
-
i
*
18
π
,
where
(6)
i
*
=
2
,
i
is
odd
,
i
*
=
1
,
i
is
even
.The layout scheme of lateral pulse jets.
(a)
Odd-numbered rings
(b)
Even-numbered ringsFigure 3
The ring frames expansion of lateral pulse jets.For the situation where all pulse jets are fired at the same time, the total force and moment are given by(7)
F
x
b
F
y
b
F
z
b
=
0
∑
i
=
1
i
=
10
‍
∑
j
=
1
j
=
18
‍
F
y
b
i
,
j
∑
i
=
1
i
=
10
‍
∑
j
=
1
j
=
18
‍
F
z
b
i
,
j
M
x
b
M
y
b
M
z
b
=
0
-
∑
i
=
1
i
=
10
‍
∑
j
=
1
j
=
18
‍
F
z
b
i
,
j
l
i
∑
i
=
1
i
=
10
‍
∑
j
=
1
j
=
18
‍
F
y
b
i
,
j
l
i
.In order to avoid the coupling between the pitch moment and yaw moment, each ring is divided into four control regions: positive pitch, negative pitch, positive yaw, and negative yaw control region, as shown in Figure4.Schematic of control regions.
(a)
Odd-numbered rings
(b)
Even-numbered ringsThe autopilot design of PAC-3 is more complicated than that of other conventional missiles which are controlled only by aerodynamic surfaces, due to the hybrid property of control inputs and the on-off property of pulse jet (the pulse jet can be fired only one time). To deal with this problem, the work [7] proposed two-step design procedures: (1) in the first step, neglect the hybrid property (or on-off property) and design the expected force and moment signals, (2) in the other step, the fire logic is derived by solving the retaliation problem of these signals, with taking into account the hybrid property (or on-off property). Different from the above traditional procedures, a novel procedure will be presented in what follows, where only one step is included. In practical applications, only a small number of jets (in a certain of rings) are activated over a finite time interval. To make the idea of the following development clear, we here consider a simple but representative situation, where no more than two rings are allowed to be fired simultaneously and no more than two jets are activated in each fired ring. Meanwhile, it should be ensured that only odd rings or even rings are fired, and the jets are fired symmetrically about the corresponding symmetry axis of each control region.Take the positive pitch control region as an example. The forces provided by the jets(
i
,
1
), (
i
,
2
), (
i
,
3
), (
i
,
17
), and (
i
,
18
) in an odd ring are given by
(8)
F
o
=
F
m
F
m
cos
⁡
π
9
F
m
cos
⁡
2
π
9
F
m
cos
⁡
2
π
9
F
m
cos
⁡
π
9
T
.Similarly, the forces provided by the jets(
i
,
1
), (
i
,
2
), (
i
,
17
), and (
i
,
18
) in an even ring are given by
(9)
F
e
=
F
m
cos
⁡
π
18
F
m
cos
⁡
π
6
F
m
cos
⁡
π
6
F
m
cos
⁡
π
18
T
,
where F
o and F
e denote the forces associated with an odd ring and an even ring, respectively.UseF
i to denote the lateral force generated by ring i. Clearly, F
i should satisfy the condition
(10)
F
1
=
F
9
,
F
3
=
F
7
,
F
2
=
F
10
,
F
4
=
F
8
.When odd ring is fired, the lateral force and moment are given by(11)
F
y
b
∈
F
m
,
2
F
m
cos
⁡
π
9
,
2
F
m
cos
⁡
2
π
9
,
2
F
m
,
4
F
m
cos
⁡
π
9
,
4
F
m
cos
⁡
2
π
9
,
M
z
b
∈
F
m
l
5
,
2
F
m
l
5
cos
⁡
π
9
,
2
F
m
l
5
cos
⁡
2
π
9
,
2
F
m
l
5
,
4
F
m
l
5
cos
⁡
π
9
,
4
F
m
l
5
cos
⁡
2
π
9
.While even ring is fired, the lateral force and moment are given by(12)
F
y
b
∈
2
F
m
cos
⁡
π
18
,
2
F
m
cos
⁡
π
6
,
4
F
m
cos
⁡
π
18
,
4
F
m
cos
⁡
π
6
M
z
b
∈
2
F
m
l
6
cos
⁡
π
18
,
2
F
m
l
6
cos
⁡
π
6
,
4
F
m
l
6
cos
⁡
π
18
,
4
F
m
l
6
cos
⁡
π
6
.NotingΔ
l is quite small, we consider l
5
≃
l
6
=
l. In order to ensure that the jets fire efficiency, the jets resulting in small moment components along z
b axis are not activated. As a result, the sets of forces and moments are, respectively, given by
(13)
U
F
y
+
=
F
m
,
2
F
m
cos
⁡
π
9
,
2
F
m
cos
⁡
π
18
,
2
F
m
,
4
F
m
cos
⁡
π
9
,
4
F
m
cos
⁡
π
18
,
(14)
U
M
y
+
=
F
m
l
,
2
F
m
l
cos
⁡
π
9
,
2
F
m
l
cos
⁡
π
18
,
2
F
m
l
,
4
F
m
l
cos
⁡
π
9
,
4
F
m
l
cos
⁡
π
18
,
where F
y
b
∈
U
F
y
+ and M
z
b
∈
U
M
y
+. During each control period, a control moment belonging to U
M
y
+ will be used as the input.The mutual interference between high-speed jet stream and air leads to lateral jet interference effect. In order to take into account this interference, force and moment amplification factorsK
F
y, K
F
z, K
M
y, and K
M
z are introduced as in [3, 4]. Then, the resulting lateral forces and moments are
(15)
F
x
b
a
F
y
b
a
F
z
b
a
=
0
F
y
b
+
K
F
y
F
y
b
F
z
b
+
K
F
z
F
z
b
,
M
x
b
a
M
y
b
a
M
z
b
a
=
0
M
y
b
+
K
M
y
M
y
b
M
z
b
+
K
M
z
M
z
b
.Remark 1.
In fact, each pulse jet can be fired only once, so the location of the fired jet cannot provide force anymore. Based on this precondition, elements of setsU
F
y
+ and U
M
y
+ will be less and less over time. In this paper, quantity change of sets’ elements is not considered to simplify the problem.


## 2.3. Attitude Control Model

Some transformation and simplification are applied to the missile model for control design. It is assumed that the missile’s mass is of a constant value. Note that the goal is to establish the angle of attack and sideslip angle. The attitude control model is given by(16)
α
˙
=
r
+
q
sin
α
tan
β
-
Q
S
C
y
α
α
+
C
y
δ
z
δ
z
cos
⁡
α
m
V
cos
⁡
β
-
F
y
b
+
K
F
y
F
y
b
cos
⁡
α
m
V
cos
⁡
β
-
G
y
cos
⁡
α
m
V
cos
⁡
β
,
β
˙
=
q
cos
⁡
α
+
(
Q
S
(
C
z
β
β
+
C
z
δ
y
δ
y
)
+
F
z
b
+
K
F
z
F
z
b
)
cos
⁡
β
m
V
+
(
Q
S
(
C
y
α
α
+
C
y
δ
z
δ
z
)
+
F
y
b
+
K
F
y
F
y
b
)
sin
α
sin
β
m
V
+
G
z
cos
⁡
β
m
V
+
G
y
sin
α
sin
β
m
V
,
q
˙
=
M
y
b
J
y
+
K
M
y
M
y
b
J
y
+
Q
S
L
m
y
β
β
J
y
+
Q
S
L
m
y
δ
y
δ
y
J
y
+
Q
S
L
m
y
q
q
J
y
,
r
˙
=
M
z
b
J
z
+
K
M
z
M
z
b
J
z
+
Q
S
L
m
z
α
α
J
z
+
Q
S
L
m
z
δ
z
δ
z
J
z
+
Q
S
L
m
z
r
r
J
z
,
where m
y
β, m
y
δ
z, m
y
q, m
z
α, m
z
δ
z, and m
z
r are aerodynamic parameters.


## 3. Mixed Logical Dynamical Model of Blended Missile

### 3.1. Piecewise Affine Model of Blended Missile

To simplify analysis, the gravity term and the channel coupling term are ignored. With (16), the missile attitude control model of pitch channel is
(17)
α
˙
=
r
-
Q
S
C
y
α
α
+
C
y
δ
z
δ
z
cos
⁡
α
m
V
-
1
+
K
F
y
F
y
b
cos
⁡
α
m
V
,
r
˙
=
(
1
+
K
M
z
)
F
y
b
l
J
z
+
Q
S
L
m
z
α
α
J
z
+
Q
S
L
m
z
δ
z
δ
z
J
z
+
Q
S
L
m
z
r
r
J
z
.Choosex
=
[
α
r
]
T as system state and u
=
[
δ
z
F
y
b
]
T as control input. The considered output is y
=
α. Then, (17) can be rewritten into the following state space form:
(18)
x
˙
=
f
x
+
g
x
u
,
y
=
1
0
x
,
where
(19)
f
(
x
)
=
r
-
Q
S
C
y
α
α
cos
⁡
α
m
V
Q
S
L
(
m
z
α
α
+
m
z
r
r
)
J
z
,
g
x
=
-
Q
S
C
y
δ
z
cos
⁡
α
m
V
-
1
+
K
F
y
cos
⁡
α
m
V
Q
S
L
m
z
δ
z
J
z
1
+
K
M
z
l
J
z
.The missile parameters are presented in Table1.Table 1
Missile’s overall parameters.
Missile’s takeoff mass [kg]
255
Missile’s full-length [m]
4.86
Missile’s diameter [m]
0.317
Distance between missile head and center of mass [m]
2.569
Moment of inertiaJ
x
,
J
y
,
J
z [kg⋅m2]
3.03, 306.3, 306.3
Lateral force provided by individual jets [N]
2200
Distance between jet ring center and center of mass [m]
1.26
Range of angle of attack [rad]
-
π
6
≤
α
≤
π
6
Range of pitch angular velocity [rad/s]
-
5
π
3
≤
r
≤
5
π
3
Range of elevator deflection [rad]
-
π
6
≤
δ
z
≤
π
6Generally, the aerodynamic coefficientsC
y
α, C
y
δ
z, m
z
α, and m
z
δ
z and the amplification factors K
F
y, K
M
z are mainly affected by the flight velocity V and the angle of attack α [11, 12]. Since the terminal guidance phase is considered in this paper, the flight time is quite short and the flight velocity of the missile can be treated as a constant. Thus, the aerodynamic coefficients and the amplification factors are mainly affected by the angle of attack α. The relationships between them are shown in Figure 5.Aerodynamic parameters as functions of angle of attack.
(a)
Relation betweenC
y
α and α
(b)
Relation betweenC
y
δ
z and α
(c)
Relation betweenm
z
α and α
(d)
Relation betweenm
z
δ
z and α
(e)
Relation betweenK
F
y and α
(f)
Relation betweenK
M
z and αIn practical application, sinceα is the main factor that leads to system nonlinearities, the system model is usually linearized if α varies in small range. As seen from Figure 5, curves of the relation between aerodynamic parameters and angle of attack can be expressed by six line segments approximately. Here, we choose α
=
-
0.37
rad, -
0.153
rad, 0,0.153
rad, and 0.37
rad as the operation points and divide the whole operation region into six subregions. As a result, the original model (17) can be converted to the following piecewise affine models:
(20)
x
˙
=
a
1
x
+
b
1
u
+
e
1
1
0
x
≤
-
0.37
a
2
x
+
b
2
u
+
e
2
-
0.37
<
1
0
x
≤
-
0.153
a
3
x
+
b
3
u
+
e
3
-
0.153
<
1
0
x
≤
0
a
4
x
+
b
4
u
+
e
4
0
<
1
0
x
≤
0.153
a
5
x
+
b
5
u
+
e
5
0.153
<
1
0
x
≤
0.37
a
6
x
+
b
6
u
+
e
6
1
0
x
>
0.37
,
(21)
y
=
1
0
x
,
where
(22)
a
i
=
∂
f
(
x
)
∂
x
x
=
x
i
0
=
a
i
11
a
i
12
a
i
21
a
i
22
,
b
i
=
g
x
i
0
=
b
i
11
b
i
12
b
i
21
b
i
22
,
e
i
=
e
i
1
e
i
2
,
a
i
11
=
∂
f
1
x
∂
α
x
=
x
i
0
=
-
Q
S
m
V
∂
C
y
α
∂
α
α
=
α
i
0
α
i
0
cos
⁡
α
i
0
+
C
y
α
(
α
i
0
)
cos
⁡
α
i
0
∂
C
y
α
∂
α
α
=
α
i
0
-
C
y
α
α
i
0
α
i
0
sin
α
i
0
,
a
i
12
=
∂
f
1
(
x
)
∂
r
x
=
x
i
0
=
1
,
a
i
21
=
∂
f
2
x
∂
α
x
=
x
i
0
=
Q
S
L
J
z
∂
m
z
α
∂
α
α
=
α
i
0
α
i
0
+
m
z
α
α
i
0
+
∂
m
z
r
∂
α
α
=
α
i
0
r
i
0
,
a
i
22
=
∂
f
2
(
x
)
∂
r
x
=
x
i
0
=
Q
S
L
J
z
m
z
α
α
i
0
,
b
i
11
=
-
Q
S
C
y
δ
z
α
i
0
cos
⁡
α
i
0
m
V
,
b
i
12
=
-
1
+
K
F
y
α
i
0
cos
⁡
α
i
0
m
V
,
b
i
21
=
Q
S
L
m
z
δ
z
(
α
i
0
)
J
z
,
b
i
22
=
(
1
+
K
M
z
(
α
i
0
)
)
l
J
z
,
e
i
1
=
a
i
11
α
i
0
+
a
i
12
r
i
0
,
e
i
2
=
a
i
21
α
i
0
+
a
i
22
r
i
0
,
where i (i
=
1,2
,
…
,
6) is the label corresponding to the ith region.From Table1 and Figure 5, we get the set of aerodynamic parameters at the point (
H
,
V
)
=
(
20
km
,
1000
m
/
s
) as shown in Table 2.Table 2
Aerodynamic parameters of pitch channel.
Aerodynamic parameter
Q
S
C
y
α
/
m
V
Q
S
C
y
δ
z
/
m
V
Q
S
L
m
z
α
/
J
z
Q
S
L
m
z
r
/
J
z
Q
S
L
m
z
δ
z
/
J
z
α
=
-
0.53
-
1.175
0.076
9.47
-
0.32
-
84.89
α
=
-
0.37
-
0.944
0.082
9.73
-
0.32
-
89.24
α
=
-
0.153
-
0.23
0.082
10.43
-
0.32
-
89.38
α
=
0
0.36
0.09
11.01
-
0.32
-
95.10
α
=
0.153
0.94
0.082
11.58
-
0.32
-
89.34
α
=
0.37
1.63
0.082
12.24
-
0.32
-
89.18Choose the sampling periodT
s
=
0.025
s. The discrete state-space expression is then given by
(23)
x
k
+
1
=
a
~
1
x
(
k
)
+
b
~
1
u
(
k
)
+
e
~
1
[
1
0
]
x
(
k
)
≤
-
0.37
a
~
2
x
(
k
)
+
b
~
2
u
(
k
)
+
e
~
2
-
0.37
<
[
1
0
]
x
(
k
)
≤
-
0.153
a
~
3
x
(
k
)
+
b
~
3
u
(
k
)
+
e
~
3
-
0.153
<
[
1
0
]
x
(
k
)
≤
0
a
~
4
x
(
k
)
+
b
~
4
u
(
k
)
+
e
~
4
0
<
[
1
0
]
x
(
k
)
≤
0.153
a
~
5
x
(
k
)
+
b
~
5
u
(
k
)
+
e
~
5
0.153
<
[
1
0
]
x
(
k
)
≤
0.37
a
~
6
x
(
k
)
+
b
~
6
u
(
k
)
+
e
~
6
[
1
0
]
x
(
k
)
>
0.37
,
(24)
y
k
=
1
0
x
k
,
where
(25)
a
~
1
=
1.04
0.025
0.22
0.995
,
a
~
2
=
1.051
0.025
0.218
0.995
,
a
~
3
=
1.023
0.0252
0.248
0.9951
,
a
~
4
=
0.9954
0.0248
0.2732
0.9954
,
a
~
5
=
0.9697
0.0245
0.2913
0.9956
,
a
~
6
=
0.9594
0.0244
0.3137
0.9959
,
e
~
1
=
-
0.0183
-
0.1123
,
e
~
2
=
-
0.0175
-
0.0792
,
e
~
3
=
-
0.0131
-
0.0376
,
e
~
4
=
0
0
,
e
~
5
=
-
0.0053
0.0454
,
e
~
6
=
-
0.0169
0.1193
,
b
~
1
=
-
0.028
1.59
×
1
0
-
6
-
2.166
1.42
×
1
0
-
4
,
b
~
2
=
-
0.03
2.18
×
1
0
-
6
-
2.224
1.87
×
1
0
-
4
,
b
~
3
=
-
0.03
2.68
×
1
0
-
6
-
2.228
2.23
×
1
0
-
4
,
b
~
4
=
-
0.032
2.77
×
1
0
-
6
-
2.371
2.32
×
1
0
-
4
,
b
~
5
=
-
0.03
2.27
×
1
0
-
6
-
2.228
1.98
×
1
0
-
4
,
b
~
6
=
-
0.029
1.91
×
1
0
-
6
-
2.224
1.73
×
1
0
-
4
.


### 3.2. Constraints Analysis

Due to the symmetry of jet configuration, the set of possible negative pitch control force is given by(26)
U
F
y
-
=
-
F
m
,
-
2
F
m
cos
⁡
π
9
,
-
2
F
m
cos
⁡
π
18
,
-
2
F
m
,
-
4
F
m
cos
⁡
π
9
,
-
4
F
m
cos
⁡
π
18
.By combining (13) and (26), we obtain the set of all possible pitch control force:
(27)
U
F
y
=
F
m
,
2
F
m
cos
⁡
π
9
,
2
F
m
cos
⁡
π
18
,
2
F
m
,
4
F
m
cos
⁡
π
9
,
4
F
m
cos
⁡
π
18
,
-
F
m
,
-
2
F
m
cos
⁡
π
9
,
-
2
F
m
cos
⁡
π
18
,
-
2
F
m
,
-
4
F
m
cos
⁡
π
9
,
-
4
F
m
cos
⁡
π
18
.Substituting the jet parameters shown in Table1 into the foregoing set yields
(28)
U
F
y
=
2200,4135,4333,4400,8269,8666
,
-
2200
,
-
4135
,
-
4333
,
-
4400
,
-
8269
,
-
8666
.Noting that the lateral forces are discrete variable, we will utilize the linear combination of logical variables to describe the lateral force. In terms of piecewise affine model (20), we introduce logical variables δ
F
i
∈
{
0,1
}, i
=
1,2
,
…
,
12 to express the lateral thrust F
y
b:
(29)
F
y
b
=
2200
δ
F
1
+
4135
δ
F
2
+
4333
δ
F
3
+
4400
δ
F
4
+
8269
δ
F
5
+
8666
δ
F
6
-
2200
δ
F
7
-
4135
δ
F
8
-
4333
δ
F
9
-
4400
δ
F
10
-
8269
δ
F
11
-
8666
δ
F
12
.In (29), the logical variables should satisfy the constraints:
(30)
∑
i
=
1
12
‍
δ
F
i
=
0
or
1
,
where 0 means that no lateral force is generated, while 1 means that the applied lateral force equals to one element of set U
F
y.Denoteu
1
=
δ
z; the control input u in model (23) can be rewritten as
(31)
u
=
u
1
F
y
b
T
.As shown in Table1, the constraints on system states and control input are
(32)
x
min
⁡
≤
x
k
≤
x
max
⁡
,
u
1
min
⁡
≤
u
1
k
≤
u
1
max
⁡
,
where x
min
⁡
=
[
-
0.53
-
5.22
]
T, x
max
⁡
=
[
0.53
5.22
]
T, u
1
min
⁡
=
-
0.53, u
1
max
⁡
=
0.53.According to (30), we have
(33)
∑
i
=
1
12
‍
δ
F
i
≤
1
which shows the constraint on the logical control input.Use logical variablesδ
i
(
k
)
∈
{
0,1
}, i
=
1,2
,
…
,
5 to describe the operation points which satisfy the constraints:
(34)
[
1
0
]
x
(
k
)
+
0.37
≤
0
⟺
δ
1
k
=
1
,
[
1
0
]
x
(
k
)
+
0.153
≤
0
⟺
δ
2
k
=
1
,
[
1
0
]
x
(
k
)
≤
0
⟺
δ
3
k
=
1
,
[
1
0
]
x
(
k
)
-
0.153
≤
0
⟺
δ
4
k
=
1
,
[
1
0
]
x
(
k
)
-
0.37
≤
0
⟺
δ
5
k
=
1
.Equation (34) can be transformed into the following equivalent mixed logical inequalities [13]:
(35)
1
0
x
k
+
0.37
≥
ɛ
+
m
1
-
ɛ
δ
1
k
,
1
0
x
k
+
0.37
≤
M
1
1
-
δ
1
k
,
1
0
x
k
+
0.153
≥
ɛ
+
m
2
-
ɛ
δ
2
k
,
1
0
x
k
+
0.153
≤
M
2
1
-
δ
2
k
,
1
0
x
k
≥
ɛ
+
m
3
-
ɛ
δ
3
k
,
1
0
x
k
≤
M
3
1
-
δ
3
k
,
1
0
x
k
-
0.153
≥
ɛ
+
m
4
-
ɛ
δ
4
k
,
1
0
x
k
-
0.153
≤
M
4
1
-
δ
4
k
,
1
0
x
k
-
0.37
≥
ɛ
+
m
5
-
ɛ
δ
5
k
,
1
0
x
k
-
0.37
≤
M
5
1
-
δ
5
k
,
where m
1
=
-
0.16, M
1
=
0.90, m
2
=
-
0.377, M
2
=
0.683, m
3
=
-
0.53, M
3
=
0.53, m
4
=
-
0.683, M
4
=
0.377, m
5
=
-
0.90, M
5
=
0.16, and ɛ
=
1
0
-
6.In addition, we introduce the auxiliary logical variablesδ
i
k
∈
0,1,  i
=
6
,
…
,
9, as follows:
(36)
δ
6
k
=
1
-
δ
1
k
δ
2
k
,
δ
7
k
=
1
-
δ
2
k
δ
3
k
,
δ
8
k
=
1
-
δ
3
k
δ
4
k
,
δ
9
k
=
1
-
δ
4
k
δ
5
k
.Withδ
1, δ
6, δ
7, δ
8, δ
9, and 1
-
δ
5, the six regions can be presented. Similarly, (36) can be transformed into the equivalent mixed logical inequalities:
(37)
-
δ
1
(
k
)
+
δ
2
(
k
)
-
δ
6
(
k
)
≤
0
,
δ
1
(
k
)
+
δ
6
(
k
)
≤
1
,
-
δ
2
(
k
)
+
δ
6
(
k
)
≤
0
,
-
δ
2
(
k
)
+
δ
3
(
k
)
-
δ
7
(
k
)
≤
0
,
δ
2
(
k
)
+
δ
7
(
k
)
≤
1
,
-
δ
3
(
k
)
+
δ
7
(
k
)
≤
0
,
-
δ
3
(
k
)
+
δ
4
(
k
)
-
δ
8
(
k
)
≤
0
,
δ
3
(
k
)
+
δ
8
(
k
)
≤
1
,
-
δ
4
(
k
)
+
δ
8
(
k
)
≤
0
,
-
δ
4
(
k
)
+
δ
5
(
k
)
-
δ
9
(
k
)
≤
0
,
δ
4
(
k
)
+
δ
9
(
k
)
≤
1
,
-
δ
5
(
k
)
+
δ
9
(
k
)
≤
0
.To describe the state space model of each region, we introduce the following auxiliary continuous variables:(38)
z
1
k
=
a
~
1
k
x
k
+
b
~
1
k
u
k
+
e
~
1
δ
1
k
,
z
2
k
=
a
~
2
k
x
k
+
b
~
2
k
u
k
+
e
~
2
δ
6
k
,
z
3
k
=
a
~
3
k
x
k
+
b
~
3
k
u
k
+
e
~
3
δ
7
k
,
z
4
k
=
a
~
4
k
x
k
+
b
~
4
k
u
k
+
e
~
4
δ
8
k
,
z
5
k
=
a
~
5
k
x
k
+
b
~
5
k
u
k
+
e
~
5
δ
9
k
,
z
6
k
=
a
~
6
k
x
k
+
b
~
6
k
u
k
+
e
~
6
1
-
δ
5
k
.Equation (38) can be converted to the equivalent mixed logical inequalities:
(39)
z
1
k
≥
a
~
1
k
x
k
+
b
~
1
k
u
k
+
e
~
1
-
M
f
1
1
-
δ
1
k
,
z
1
k
≤
a
~
1
k
x
k
+
b
~
1
k
u
k
+
e
~
1
-
m
f
1
1
-
δ
1
k
,
z
1
k
≥
m
f
1
δ
1
k
,
z
1
k
≤
M
f
1
δ
1
k
,
z
2
k
≥
a
~
2
k
x
k
+
b
~
2
k
u
k
+
e
~
2
-
M
f
2
1
-
δ
6
k
,
z
2
k
≤
a
~
2
k
x
k
+
b
~
2
k
u
k
+
e
~
2
-
m
f
2
1
-
δ
6
k
,
z
2
k
≥
m
f
2
δ
6
k
,
z
2
k
≤
M
f
2
δ
6
k
,
z
3
k
≥
a
~
3
k
x
k
+
b
~
3
k
u
k
+
e
~
3
-
M
f
3
1
-
δ
7
k
,
z
3
k
≤
a
~
3
k
x
k
+
b
~
3
k
u
k
+
e
~
3
-
m
f
3
1
-
δ
7
k
,
z
3
k
≥
m
f
3
δ
7
k
,
z
3
k
≤
M
f
3
δ
7
k
,
z
4
k
≥
a
~
4
k
x
k
+
b
~
4
k
u
k
+
e
~
4
-
M
f
4
1
-
δ
8
k
,
z
4
k
≤
a
~
4
k
x
k
+
b
~
4
k
u
k
+
e
~
4
-
m
f
4
1
-
δ
8
k
,
z
4
k
≥
m
f
4
δ
8
k
,
z
4
k
≤
M
f
4
δ
8
k
,
z
5
k
≥
a
~
5
k
x
k
+
b
~
5
k
u
k
+
e
~
5
-
M
f
5
1
-
δ
9
k
,
z
5
k
≤
a
~
5
k
x
k
+
b
~
5
k
u
k
+
e
~
5
-
m
f
5
1
-
δ
9
k
,
z
5
k
≥
m
f
5
δ
9
k
,
z
5
k
≤
M
f
5
δ
9
k
,
z
6
k
≥
a
~
6
k
x
k
+
b
~
6
k
u
k
+
e
~
6
-
M
f
5
δ
5
k
,
z
6
k
≤
a
~
6
k
x
k
+
b
~
6
k
u
k
+
e
~
6
-
m
f
5
δ
5
k
,
z
6
k
≥
m
f
6
1
-
δ
5
k
,
z
6
k
≤
M
f
6
1
-
δ
5
k
,
where M
f
1
=
[
0.73
10.84
]
T, m
f
1
=
[
-
0.77
-
11.06
]
T, M
f
2
=
[
0.76
11.66
]
T, m
f
2
=
[
-
0.79
-
12.56
]
T, M
f
3
=
[
0.77
12.73
]
T, m
f
3
=
[
-
0.78
-
12.80
]
T, M
f
4
=
[
0.76
14.03
]
T, m
f
4
=
[
-
0.76
-
14.03
]
T, M
f
5
=
[
0.725
12.91
]
T, m
f
5
=
[
-
0.736
-
12.05
]
T, M
f
6
=
[
0.696
11.54
]
T, and m
f
6
=
[
-
0.73
-
11.30
]
T.Then, the whole MLD model of the missile is given by(40)
x
(
k
+
1
)
=
∑
i
=
1
6
‍
z
i
(
k
)
y
(
k
)
=
[
1
0
]
x
(
k
)
s
.
t
.
28
,
30
,
32
,
34
,
36
,
38
,
where (29), (31), (33), (35), (37), and (39) represent all the constraint inequalities.


## 3.1. Piecewise Affine Model of Blended Missile

To simplify analysis, the gravity term and the channel coupling term are ignored. With (16), the missile attitude control model of pitch channel is
(17)
α
˙
=
r
-
Q
S
C
y
α
α
+
C
y
δ
z
δ
z
cos
⁡
α
m
V
-
1
+
K
F
y
F
y
b
cos
⁡
α
m
V
,
r
˙
=
(
1
+
K
M
z
)
F
y
b
l
J
z
+
Q
S
L
m
z
α
α
J
z
+
Q
S
L
m
z
δ
z
δ
z
J
z
+
Q
S
L
m
z
r
r
J
z
.Choosex
=
[
α
r
]
T as system state and u
=
[
δ
z
F
y
b
]
T as control input. The considered output is y
=
α. Then, (17) can be rewritten into the following state space form:
(18)
x
˙
=
f
x
+
g
x
u
,
y
=
1
0
x
,
where
(19)
f
(
x
)
=
r
-
Q
S
C
y
α
α
cos
⁡
α
m
V
Q
S
L
(
m
z
α
α
+
m
z
r
r
)
J
z
,
g
x
=
-
Q
S
C
y
δ
z
cos
⁡
α
m
V
-
1
+
K
F
y
cos
⁡
α
m
V
Q
S
L
m
z
δ
z
J
z
1
+
K
M
z
l
J
z
.The missile parameters are presented in Table1.Table 1
Missile’s overall parameters.
Missile’s takeoff mass [kg]
255
Missile’s full-length [m]
4.86
Missile’s diameter [m]
0.317
Distance between missile head and center of mass [m]
2.569
Moment of inertiaJ
x
,
J
y
,
J
z [kg⋅m2]
3.03, 306.3, 306.3
Lateral force provided by individual jets [N]
2200
Distance between jet ring center and center of mass [m]
1.26
Range of angle of attack [rad]
-
π
6
≤
α
≤
π
6
Range of pitch angular velocity [rad/s]
-
5
π
3
≤
r
≤
5
π
3
Range of elevator deflection [rad]
-
π
6
≤
δ
z
≤
π
6Generally, the aerodynamic coefficientsC
y
α, C
y
δ
z, m
z
α, and m
z
δ
z and the amplification factors K
F
y, K
M
z are mainly affected by the flight velocity V and the angle of attack α [11, 12]. Since the terminal guidance phase is considered in this paper, the flight time is quite short and the flight velocity of the missile can be treated as a constant. Thus, the aerodynamic coefficients and the amplification factors are mainly affected by the angle of attack α. The relationships between them are shown in Figure 5.Aerodynamic parameters as functions of angle of attack.
(a)
Relation betweenC
y
α and α
(b)
Relation betweenC
y
δ
z and α
(c)
Relation betweenm
z
α and α
(d)
Relation betweenm
z
δ
z and α
(e)
Relation betweenK
F
y and α
(f)
Relation betweenK
M
z and αIn practical application, sinceα is the main factor that leads to system nonlinearities, the system model is usually linearized if α varies in small range. As seen from Figure 5, curves of the relation between aerodynamic parameters and angle of attack can be expressed by six line segments approximately. Here, we choose α
=
-
0.37
rad, -
0.153
rad, 0,0.153
rad, and 0.37
rad as the operation points and divide the whole operation region into six subregions. As a result, the original model (17) can be converted to the following piecewise affine models:
(20)
x
˙
=
a
1
x
+
b
1
u
+
e
1
1
0
x
≤
-
0.37
a
2
x
+
b
2
u
+
e
2
-
0.37
<
1
0
x
≤
-
0.153
a
3
x
+
b
3
u
+
e
3
-
0.153
<
1
0
x
≤
0
a
4
x
+
b
4
u
+
e
4
0
<
1
0
x
≤
0.153
a
5
x
+
b
5
u
+
e
5
0.153
<
1
0
x
≤
0.37
a
6
x
+
b
6
u
+
e
6
1
0
x
>
0.37
,
(21)
y
=
1
0
x
,
where
(22)
a
i
=
∂
f
(
x
)
∂
x
x
=
x
i
0
=
a
i
11
a
i
12
a
i
21
a
i
22
,
b
i
=
g
x
i
0
=
b
i
11
b
i
12
b
i
21
b
i
22
,
e
i
=
e
i
1
e
i
2
,
a
i
11
=
∂
f
1
x
∂
α
x
=
x
i
0
=
-
Q
S
m
V
∂
C
y
α
∂
α
α
=
α
i
0
α
i
0
cos
⁡
α
i
0
+
C
y
α
(
α
i
0
)
cos
⁡
α
i
0
∂
C
y
α
∂
α
α
=
α
i
0
-
C
y
α
α
i
0
α
i
0
sin
α
i
0
,
a
i
12
=
∂
f
1
(
x
)
∂
r
x
=
x
i
0
=
1
,
a
i
21
=
∂
f
2
x
∂
α
x
=
x
i
0
=
Q
S
L
J
z
∂
m
z
α
∂
α
α
=
α
i
0
α
i
0
+
m
z
α
α
i
0
+
∂
m
z
r
∂
α
α
=
α
i
0
r
i
0
,
a
i
22
=
∂
f
2
(
x
)
∂
r
x
=
x
i
0
=
Q
S
L
J
z
m
z
α
α
i
0
,
b
i
11
=
-
Q
S
C
y
δ
z
α
i
0
cos
⁡
α
i
0
m
V
,
b
i
12
=
-
1
+
K
F
y
α
i
0
cos
⁡
α
i
0
m
V
,
b
i
21
=
Q
S
L
m
z
δ
z
(
α
i
0
)
J
z
,
b
i
22
=
(
1
+
K
M
z
(
α
i
0
)
)
l
J
z
,
e
i
1
=
a
i
11
α
i
0
+
a
i
12
r
i
0
,
e
i
2
=
a
i
21
α
i
0
+
a
i
22
r
i
0
,
where i (i
=
1,2
,
…
,
6) is the label corresponding to the ith region.From Table1 and Figure 5, we get the set of aerodynamic parameters at the point (
H
,
V
)
=
(
20
km
,
1000
m
/
s
) as shown in Table 2.Table 2
Aerodynamic parameters of pitch channel.
Aerodynamic parameter
Q
S
C
y
α
/
m
V
Q
S
C
y
δ
z
/
m
V
Q
S
L
m
z
α
/
J
z
Q
S
L
m
z
r
/
J
z
Q
S
L
m
z
δ
z
/
J
z
α
=
-
0.53
-
1.175
0.076
9.47
-
0.32
-
84.89
α
=
-
0.37
-
0.944
0.082
9.73
-
0.32
-
89.24
α
=
-
0.153
-
0.23
0.082
10.43
-
0.32
-
89.38
α
=
0
0.36
0.09
11.01
-
0.32
-
95.10
α
=
0.153
0.94
0.082
11.58
-
0.32
-
89.34
α
=
0.37
1.63
0.082
12.24
-
0.32
-
89.18Choose the sampling periodT
s
=
0.025
s. The discrete state-space expression is then given by
(23)
x
k
+
1
=
a
~
1
x
(
k
)
+
b
~
1
u
(
k
)
+
e
~
1
[
1
0
]
x
(
k
)
≤
-
0.37
a
~
2
x
(
k
)
+
b
~
2
u
(
k
)
+
e
~
2
-
0.37
<
[
1
0
]
x
(
k
)
≤
-
0.153
a
~
3
x
(
k
)
+
b
~
3
u
(
k
)
+
e
~
3
-
0.153
<
[
1
0
]
x
(
k
)
≤
0
a
~
4
x
(
k
)
+
b
~
4
u
(
k
)
+
e
~
4
0
<
[
1
0
]
x
(
k
)
≤
0.153
a
~
5
x
(
k
)
+
b
~
5
u
(
k
)
+
e
~
5
0.153
<
[
1
0
]
x
(
k
)
≤
0.37
a
~
6
x
(
k
)
+
b
~
6
u
(
k
)
+
e
~
6
[
1
0
]
x
(
k
)
>
0.37
,
(24)
y
k
=
1
0
x
k
,
where
(25)
a
~
1
=
1.04
0.025
0.22
0.995
,
a
~
2
=
1.051
0.025
0.218
0.995
,
a
~
3
=
1.023
0.0252
0.248
0.9951
,
a
~
4
=
0.9954
0.0248
0.2732
0.9954
,
a
~
5
=
0.9697
0.0245
0.2913
0.9956
,
a
~
6
=
0.9594
0.0244
0.3137
0.9959
,
e
~
1
=
-
0.0183
-
0.1123
,
e
~
2
=
-
0.0175
-
0.0792
,
e
~
3
=
-
0.0131
-
0.0376
,
e
~
4
=
0
0
,
e
~
5
=
-
0.0053
0.0454
,
e
~
6
=
-
0.0169
0.1193
,
b
~
1
=
-
0.028
1.59
×
1
0
-
6
-
2.166
1.42
×
1
0
-
4
,
b
~
2
=
-
0.03
2.18
×
1
0
-
6
-
2.224
1.87
×
1
0
-
4
,
b
~
3
=
-
0.03
2.68
×
1
0
-
6
-
2.228
2.23
×
1
0
-
4
,
b
~
4
=
-
0.032
2.77
×
1
0
-
6
-
2.371
2.32
×
1
0
-
4
,
b
~
5
=
-
0.03
2.27
×
1
0
-
6
-
2.228
1.98
×
1
0
-
4
,
b
~
6
=
-
0.029
1.91
×
1
0
-
6
-
2.224
1.73
×
1
0
-
4
.


## 3.2. Constraints Analysis

Due to the symmetry of jet configuration, the set of possible negative pitch control force is given by(26)
U
F
y
-
=
-
F
m
,
-
2
F
m
cos
⁡
π
9
,
-
2
F
m
cos
⁡
π
18
,
-
2
F
m
,
-
4
F
m
cos
⁡
π
9
,
-
4
F
m
cos
⁡
π
18
.By combining (13) and (26), we obtain the set of all possible pitch control force:
(27)
U
F
y
=
F
m
,
2
F
m
cos
⁡
π
9
,
2
F
m
cos
⁡
π
18
,
2
F
m
,
4
F
m
cos
⁡
π
9
,
4
F
m
cos
⁡
π
18
,
-
F
m
,
-
2
F
m
cos
⁡
π
9
,
-
2
F
m
cos
⁡
π
18
,
-
2
F
m
,
-
4
F
m
cos
⁡
π
9
,
-
4
F
m
cos
⁡
π
18
.Substituting the jet parameters shown in Table1 into the foregoing set yields
(28)
U
F
y
=
2200,4135,4333,4400,8269,8666
,
-
2200
,
-
4135
,
-
4333
,
-
4400
,
-
8269
,
-
8666
.Noting that the lateral forces are discrete variable, we will utilize the linear combination of logical variables to describe the lateral force. In terms of piecewise affine model (20), we introduce logical variables δ
F
i
∈
{
0,1
}, i
=
1,2
,
…
,
12 to express the lateral thrust F
y
b:
(29)
F
y
b
=
2200
δ
F
1
+
4135
δ
F
2
+
4333
δ
F
3
+
4400
δ
F
4
+
8269
δ
F
5
+
8666
δ
F
6
-
2200
δ
F
7
-
4135
δ
F
8
-
4333
δ
F
9
-
4400
δ
F
10
-
8269
δ
F
11
-
8666
δ
F
12
.In (29), the logical variables should satisfy the constraints:
(30)
∑
i
=
1
12
‍
δ
F
i
=
0
or
1
,
where 0 means that no lateral force is generated, while 1 means that the applied lateral force equals to one element of set U
F
y.Denoteu
1
=
δ
z; the control input u in model (23) can be rewritten as
(31)
u
=
u
1
F
y
b
T
.As shown in Table1, the constraints on system states and control input are
(32)
x
min
⁡
≤
x
k
≤
x
max
⁡
,
u
1
min
⁡
≤
u
1
k
≤
u
1
max
⁡
,
where x
min
⁡
=
[
-
0.53
-
5.22
]
T, x
max
⁡
=
[
0.53
5.22
]
T, u
1
min
⁡
=
-
0.53, u
1
max
⁡
=
0.53.According to (30), we have
(33)
∑
i
=
1
12
‍
δ
F
i
≤
1
which shows the constraint on the logical control input.Use logical variablesδ
i
(
k
)
∈
{
0,1
}, i
=
1,2
,
…
,
5 to describe the operation points which satisfy the constraints:
(34)
[
1
0
]
x
(
k
)
+
0.37
≤
0
⟺
δ
1
k
=
1
,
[
1
0
]
x
(
k
)
+
0.153
≤
0
⟺
δ
2
k
=
1
,
[
1
0
]
x
(
k
)
≤
0
⟺
δ
3
k
=
1
,
[
1
0
]
x
(
k
)
-
0.153
≤
0
⟺
δ
4
k
=
1
,
[
1
0
]
x
(
k
)
-
0.37
≤
0
⟺
δ
5
k
=
1
.Equation (34) can be transformed into the following equivalent mixed logical inequalities [13]:
(35)
1
0
x
k
+
0.37
≥
ɛ
+
m
1
-
ɛ
δ
1
k
,
1
0
x
k
+
0.37
≤
M
1
1
-
δ
1
k
,
1
0
x
k
+
0.153
≥
ɛ
+
m
2
-
ɛ
δ
2
k
,
1
0
x
k
+
0.153
≤
M
2
1
-
δ
2
k
,
1
0
x
k
≥
ɛ
+
m
3
-
ɛ
δ
3
k
,
1
0
x
k
≤
M
3
1
-
δ
3
k
,
1
0
x
k
-
0.153
≥
ɛ
+
m
4
-
ɛ
δ
4
k
,
1
0
x
k
-
0.153
≤
M
4
1
-
δ
4
k
,
1
0
x
k
-
0.37
≥
ɛ
+
m
5
-
ɛ
δ
5
k
,
1
0
x
k
-
0.37
≤
M
5
1
-
δ
5
k
,
where m
1
=
-
0.16, M
1
=
0.90, m
2
=
-
0.377, M
2
=
0.683, m
3
=
-
0.53, M
3
=
0.53, m
4
=
-
0.683, M
4
=
0.377, m
5
=
-
0.90, M
5
=
0.16, and ɛ
=
1
0
-
6.In addition, we introduce the auxiliary logical variablesδ
i
k
∈
0,1,  i
=
6
,
…
,
9, as follows:
(36)
δ
6
k
=
1
-
δ
1
k
δ
2
k
,
δ
7
k
=
1
-
δ
2
k
δ
3
k
,
δ
8
k
=
1
-
δ
3
k
δ
4
k
,
δ
9
k
=
1
-
δ
4
k
δ
5
k
.Withδ
1, δ
6, δ
7, δ
8, δ
9, and 1
-
δ
5, the six regions can be presented. Similarly, (36) can be transformed into the equivalent mixed logical inequalities:
(37)
-
δ
1
(
k
)
+
δ
2
(
k
)
-
δ
6
(
k
)
≤
0
,
δ
1
(
k
)
+
δ
6
(
k
)
≤
1
,
-
δ
2
(
k
)
+
δ
6
(
k
)
≤
0
,
-
δ
2
(
k
)
+
δ
3
(
k
)
-
δ
7
(
k
)
≤
0
,
δ
2
(
k
)
+
δ
7
(
k
)
≤
1
,
-
δ
3
(
k
)
+
δ
7
(
k
)
≤
0
,
-
δ
3
(
k
)
+
δ
4
(
k
)
-
δ
8
(
k
)
≤
0
,
δ
3
(
k
)
+
δ
8
(
k
)
≤
1
,
-
δ
4
(
k
)
+
δ
8
(
k
)
≤
0
,
-
δ
4
(
k
)
+
δ
5
(
k
)
-
δ
9
(
k
)
≤
0
,
δ
4
(
k
)
+
δ
9
(
k
)
≤
1
,
-
δ
5
(
k
)
+
δ
9
(
k
)
≤
0
.To describe the state space model of each region, we introduce the following auxiliary continuous variables:(38)
z
1
k
=
a
~
1
k
x
k
+
b
~
1
k
u
k
+
e
~
1
δ
1
k
,
z
2
k
=
a
~
2
k
x
k
+
b
~
2
k
u
k
+
e
~
2
δ
6
k
,
z
3
k
=
a
~
3
k
x
k
+
b
~
3
k
u
k
+
e
~
3
δ
7
k
,
z
4
k
=
a
~
4
k
x
k
+
b
~
4
k
u
k
+
e
~
4
δ
8
k
,
z
5
k
=
a
~
5
k
x
k
+
b
~
5
k
u
k
+
e
~
5
δ
9
k
,
z
6
k
=
a
~
6
k
x
k
+
b
~
6
k
u
k
+
e
~
6
1
-
δ
5
k
.Equation (38) can be converted to the equivalent mixed logical inequalities:
(39)
z
1
k
≥
a
~
1
k
x
k
+
b
~
1
k
u
k
+
e
~
1
-
M
f
1
1
-
δ
1
k
,
z
1
k
≤
a
~
1
k
x
k
+
b
~
1
k
u
k
+
e
~
1
-
m
f
1
1
-
δ
1
k
,
z
1
k
≥
m
f
1
δ
1
k
,
z
1
k
≤
M
f
1
δ
1
k
,
z
2
k
≥
a
~
2
k
x
k
+
b
~
2
k
u
k
+
e
~
2
-
M
f
2
1
-
δ
6
k
,
z
2
k
≤
a
~
2
k
x
k
+
b
~
2
k
u
k
+
e
~
2
-
m
f
2
1
-
δ
6
k
,
z
2
k
≥
m
f
2
δ
6
k
,
z
2
k
≤
M
f
2
δ
6
k
,
z
3
k
≥
a
~
3
k
x
k
+
b
~
3
k
u
k
+
e
~
3
-
M
f
3
1
-
δ
7
k
,
z
3
k
≤
a
~
3
k
x
k
+
b
~
3
k
u
k
+
e
~
3
-
m
f
3
1
-
δ
7
k
,
z
3
k
≥
m
f
3
δ
7
k
,
z
3
k
≤
M
f
3
δ
7
k
,
z
4
k
≥
a
~
4
k
x
k
+
b
~
4
k
u
k
+
e
~
4
-
M
f
4
1
-
δ
8
k
,
z
4
k
≤
a
~
4
k
x
k
+
b
~
4
k
u
k
+
e
~
4
-
m
f
4
1
-
δ
8
k
,
z
4
k
≥
m
f
4
δ
8
k
,
z
4
k
≤
M
f
4
δ
8
k
,
z
5
k
≥
a
~
5
k
x
k
+
b
~
5
k
u
k
+
e
~
5
-
M
f
5
1
-
δ
9
k
,
z
5
k
≤
a
~
5
k
x
k
+
b
~
5
k
u
k
+
e
~
5
-
m
f
5
1
-
δ
9
k
,
z
5
k
≥
m
f
5
δ
9
k
,
z
5
k
≤
M
f
5
δ
9
k
,
z
6
k
≥
a
~
6
k
x
k
+
b
~
6
k
u
k
+
e
~
6
-
M
f
5
δ
5
k
,
z
6
k
≤
a
~
6
k
x
k
+
b
~
6
k
u
k
+
e
~
6
-
m
f
5
δ
5
k
,
z
6
k
≥
m
f
6
1
-
δ
5
k
,
z
6
k
≤
M
f
6
1
-
δ
5
k
,
where M
f
1
=
[
0.73
10.84
]
T, m
f
1
=
[
-
0.77
-
11.06
]
T, M
f
2
=
[
0.76
11.66
]
T, m
f
2
=
[
-
0.79
-
12.56
]
T, M
f
3
=
[
0.77
12.73
]
T, m
f
3
=
[
-
0.78
-
12.80
]
T, M
f
4
=
[
0.76
14.03
]
T, m
f
4
=
[
-
0.76
-
14.03
]
T, M
f
5
=
[
0.725
12.91
]
T, m
f
5
=
[
-
0.736
-
12.05
]
T, M
f
6
=
[
0.696
11.54
]
T, and m
f
6
=
[
-
0.73
-
11.30
]
T.Then, the whole MLD model of the missile is given by(40)
x
(
k
+
1
)
=
∑
i
=
1
6
‍
z
i
(
k
)
y
(
k
)
=
[
1
0
]
x
(
k
)
s
.
t
.
28
,
30
,
32
,
34
,
36
,
38
,
where (29), (31), (33), (35), (37), and (39) represent all the constraint inequalities.


## 4. Autopilot Design Using Hybrid MPC Method

In general, we expect that the outputy tracks its command y
c as fast as possible with a small amount of fuel consumption. Motivated by this observation, we consider the following optimization problem:
(41)
J
*
=

min
⁡
u
k
,
u
k
+
1
,
δ
k
,
δ
k
+
1
∣
k
,
z
k
,
z
k
+
1
∣
k
∑
i
=
1
N
‍
y
k
+
i
∣
k
y
c
k
+
i
u
(
k
+
i
)
R
2
&
11111111111111111111111111111
-
y
c
k
+
i
Q
2
&
1111111111111111111111111111
+
u
(
k
+
i
)
R
2
s
.
t
.
MLD
model
(
39
)
s
.
t
.
u
1
min
⁡
≤
u
1
k
,
u
1
k
+
1
≤
u
1
max
⁡
s
.
t
.
x
min
⁡
≤
x
k
,
x
k
+
1
≤
x
max
⁡
,
where y
c is the command for angle of attack, y
(
k
+
i
∣
k
) is predictive value of angle of attack, N denotes the predictive horizon, Q and R are the weighting matrices, and R
=
diag
⁡
(
R
δ
z
,
R
F
y
). R
δ
z, R
F
y are weighted coefficients of aerodynamic control surfaces and lateral pulse jets, respectively. The ratio of R
δ
z and R
F
y represents the control allocation between aerodynamic control surfaces and lateral pulse jets. When R
F
y
/
R
δ
z is increased, the requirement for lateral force will be declined; that is, the fuel consumption can be lessened. In addition, the deduction can be verified from the simulation results of Case 1 and Case 2.Different control allocation between the dual actuators (aerodynamic control surfaces and lateral pulse jets) can be obtained by setting differentQ and R. Since both logical and continuous variables are involved in (41), the above on-line optimization problem is a mixed integer quadratic programming (MIQP) problem which can be solved by using the hybrid MPC toolbox of MATLAB [14, 15].Remark 2.
In an explicit MPC controller, the main factors affecting the number of subregions include the system dimension, the predictive horizon, and the number of constrains. In this paper, the blended missile with aerodynamic control surfaces and lateral jets is investigated. To deal with the discrete property of the lateral jet forces, some logical variables are introduced for MPC controller design, which increases the number of subregions inevitably.


## 5. Numerical Simulations

### 5.1. Case 1

Numerical simulations are performed to verify the feasibility of the proposed method in this section. The involved parameters are given in Table1. The weight matrices are set as Q
=
10, R
δ
z
=
0.025, and R
F
y
=
diag
⁡
(
0.005,0.005
,
…
,
0.005
). Suppose the initial state is α
0
=
0, r
0
=
0. In order to avoid excessive computational load, we choose predictive horizon N
=
2 and control horizon M
=
2. The command for angle of attack is α
c
=
0.2
rad.Explicit form of the optimal controller is provided by hybrid toolbox, which is not given here for simplicity. Instead, we present the partition results associated with the controller in Figure6, where 4190 subregions are obtained and, for each subregion, a control law of the form u
=
K
x
+
b is given. As a benefit of using explicit MPC, the controller parameters K and b for all the subregions are obtained simultaneously once the optimization problem is solved. This property is obviously different from that of gain-scheduled controller (where the controller parameters for different regions are separately determined).Figure 6
The partition associated with the explicit controller.However, by a simple calculation, we conclude that these subregions do not require too much storage space (the needed storage space is approximately14
×
8
×
4190
=
450
KB). Moreover, it should be noted that when the algorithm is implemented in practice, much more time will be needed to identify which subregion the current states enter into. This is another important factor that affects the computational cost of the proposed algorithm.The simulation results are shown in Figures7–9. It is seen from Figure 7 that the actual angle of attack tracks the command in less than 0.1 seconds. As shown in Figures 8 and 9, the aerodynamic force converges to a constant as time increases while the jets are only activated during the beginning period (when the tracking error is obvious).Figure 7
The response of angle of attack.Figure 8
The response of elevator deflection.Figure 9
The response of lateral force.


### 5.2. Case 2

The results for the case with a different weightR
δ
z
=
0.005 are shown in Figures 10–12. It is seen that a different control allocation result is obtained.Figure 10
The response of angle of attack.Figure 11
The response of elevator deflection.Figure 12
The response of lateral force.


### 5.3. Case 3

In this case, a different commandα
c
=
0.2
cos
⁡
(
0.5
π
t
)
rad is considered and the weight matrices are the same as the ones for Case 1. The simulation results are shown in Figures 13–15.Figure 13
The response of angle of attack.Figure 14
The response of elevator deflection.Figure 15
The response of lateral force.It is seen from Figure13 that asymptotic tracking is also achieved. The method proposed in this paper can realize both fast tracking command and control allocation. It is seen from Figure 15 that the jets are activated, while the tracking error is obvious or command is varying rapidly.


## 5.1. Case 1

Numerical simulations are performed to verify the feasibility of the proposed method in this section. The involved parameters are given in Table1. The weight matrices are set as Q
=
10, R
δ
z
=
0.025, and R
F
y
=
diag
⁡
(
0.005,0.005
,
…
,
0.005
). Suppose the initial state is α
0
=
0, r
0
=
0. In order to avoid excessive computational load, we choose predictive horizon N
=
2 and control horizon M
=
2. The command for angle of attack is α
c
=
0.2
rad.Explicit form of the optimal controller is provided by hybrid toolbox, which is not given here for simplicity. Instead, we present the partition results associated with the controller in Figure6, where 4190 subregions are obtained and, for each subregion, a control law of the form u
=
K
x
+
b is given. As a benefit of using explicit MPC, the controller parameters K and b for all the subregions are obtained simultaneously once the optimization problem is solved. This property is obviously different from that of gain-scheduled controller (where the controller parameters for different regions are separately determined).Figure 6
The partition associated with the explicit controller.However, by a simple calculation, we conclude that these subregions do not require too much storage space (the needed storage space is approximately14
×
8
×
4190
=
450
KB). Moreover, it should be noted that when the algorithm is implemented in practice, much more time will be needed to identify which subregion the current states enter into. This is another important factor that affects the computational cost of the proposed algorithm.The simulation results are shown in Figures7–9. It is seen from Figure 7 that the actual angle of attack tracks the command in less than 0.1 seconds. As shown in Figures 8 and 9, the aerodynamic force converges to a constant as time increases while the jets are only activated during the beginning period (when the tracking error is obvious).Figure 7
The response of angle of attack.Figure 8
The response of elevator deflection.Figure 9
The response of lateral force.


## 5.2. Case 2

The results for the case with a different weightR
δ
z
=
0.005 are shown in Figures 10–12. It is seen that a different control allocation result is obtained.Figure 10
The response of angle of attack.Figure 11
The response of elevator deflection.Figure 12
The response of lateral force.


## 5.3. Case 3

In this case, a different commandα
c
=
0.2
cos
⁡
(
0.5
π
t
)
rad is considered and the weight matrices are the same as the ones for Case 1. The simulation results are shown in Figures 13–15.Figure 13
The response of angle of attack.Figure 14
The response of elevator deflection.Figure 15
The response of lateral force.It is seen from Figure13 that asymptotic tracking is also achieved. The method proposed in this paper can realize both fast tracking command and control allocation. It is seen from Figure 15 that the jets are activated, while the tracking error is obvious or command is varying rapidly.


## 6. Conclusion

An autopilot design method for a missile with aerodynamic control surfaces and lateral jets is presented in this paper. The nonlinear attitude control model is reduced to an MLD model. Meanwhile, the lateral force is described as linear combination of logical variables due to the discrete values of lateral force. Then the whole MLD model of attitude control system is derived. Autopilot design is accomplished using hybrid MPC method. By setting the related weighted coefficients in index function, the control allocation is obtained. Moreover, numerical simulations are performed under the different conditions, the performance of tracking the attitude command and control allocation is verified, and the explicit form of the control law can be obtained as well.



---

*Source: 718036-2015-02-08.xml*