# Velocity Control for Coning Motion Missile System Using Direct Discretization Method

Authors: Rui-sheng  Sun; Chao Ming
Journal: Discrete Dynamics in Nature and Society
 (2015)

---

## Abstract

This paper presents a new coning motion control methodology, which takes into account the terminal speed constraint to design the velocity control system for a missile. By using a direct discretization method to transform the optimal control problem into a nonlinear dynamic programming problem, the optimal trajectory and velocity profile are obtained to satisfy the design index requirement. In order to perform the velocity control, a virtual moving target is proposed for the missile to chase along the optimized trajectory. Consequently, after building velocity control model, a velocity control law and control parameters of the coning motion are completed through the dynamic inversion theory. The simulation results suggest that the proposed control law has a good performance and could be applied to the guidance for the missile with terminal speed control constraint.


---

## Body

## 1. Introduction

Modern warfare demands increasingly high performance requirements for missile, for example, reentry warhead, cluster bombs, and rocket-assisted torpedo. It needs not only high guidance precision, but also desired velocity in terminal phase. The reason for this is threefold. The first one is the conceptual requirement in terminal phase. For instance, a rocket-assisted torpedo should fly not too fast to guarantee that the matrix and the subtorpedo can separate securely. The second is the requirement of the large launch window and large envelope. For a missile system, it should not only get as large range as possible, but also hit the enough short-range targets, which requires the missile to consume superfluous kinetic energy in a timely manner. The third is to ensure the performance of control system. For example, the radar signals can fail in transfer due to the high-speed motion of radar-guided missile, which makes reentry warhead surface be surrounded by a serious aerodynamic heating of plasma. Consequently, it is necessary to control the deceleration of the missile flight speed in a proper time, which has become a focus for researchers.Generally speaking, there are three kinds of approaches to control the velocity of missile when flying in the air. Firstly, we can control the velocity by changing the thrust force magnitude. Enomoto et al. [1] introduced a velocity control system for a leader-following UAV through changing the thrust force magnitude and using the dynamic inversion with the two-time scale approach. However, this control approach can only be suitable for the missile with propulsion system still on.Secondly, according to the flight status and the terminal velocity constraint, the missile velocity can be governed by using the trajectory optimization method. By using this approach, the flight velocity is controlled by changing the flight height profile of the missile, which can change its drag force due to its changing dynamic pressure. In [2], a defined problem of optimal control was transformed into a two-point boundary problem, through applying the Pontryagin minimum principle. In this work, the optimal control and the optimal trajectory corresponding to this control at a skip three-dimensional entry of space vehicle into the planetary atmosphere were determined related to the obtaining of the maximum terminal velocity. Saraf et al. [3] presented an entry control algorithm for future space transportation vehicles through tracking the reference drag acceleration and heading angle profiles so as to satisfy all entry constraints including flight velocity control. Bruyère et al. [4, 5] designed a sideslip velocity autopilot for a model of tactical missile in order to meet the requirements of sideslip velocity over full flight envelope.Finally, the missile velocity can be controlled by means of the coning motion method. It is a type of control technique that the missile axis movement is tapered at a fixed angle of rotation around the velocity vector, which will produce a velocity-centered conical surface, a coning motion [6] that can make the missile increase induced drag in order to adjust the flight velocity through producing induced angle of attack. Obviously, the coning motion control is an effective approach to control the flight velocity of missile [7–9]. Song [10] presented a new velocity control method of coning motion through simulating the reentry to improve flight control accuracy and control robustness of the terminal speed of the reentry vehicle when antidesigning the speed control method for Pershing II. Reference [11] put forward a velocity magnitude control program which separates deceleration motion into two different forms in the capability range of warhead guidance and control systems, in order to meet the requirement of fall velocity of a reentry maneuvering warhead and to avoid designing a complicated ideal velocity curve. In the field of noncontrolled rocket, Mao et al. did some work for the coning motion on the producing mechanism [12], motion characteristics and stability analysis methods [8, 13], optimal control [14], and the approach of reducing and avoiding coning motion [15].Consequently, in order to design the velocity control system for a missile with the constraint of the terminal velocity, a deep understanding of the interaction among flight mechanics, optimizations, and control is necessary. The aim of this research is to propose a design methodology for the velocity control system by utilizing the coning motion based on a direct discretization and dynamic inversion method. The remaining of this paper is organized as follows. Section2 describes the problem formulation and design scheme for velocity control. In Section 3, a direct discretization method is deduced to develop the standard trajectory optimization and the velocity control profile, and the velocity control methodology is presented in detail through the coning motion technique. Moreover, the structure and parameter design of velocity controller are explored for a coning motion missile. In Section 4, a simulation case is demonstrated to govern the terminal velocity of the missile through using the coning motion algorithm. In order to illustrate the velocity control performance, a traditional control method is also carried out for comparative study. Finally, some conclusions are presented in Section 5.


## 2. Problem Formulation

### 2.1. Motion Equations

For the optimal flight velocity control problem under consideration, a mathematical model based on point mass dynamics can be used for determining the trajectory and velocity, yielding(1)mV˙M=-D-mgsin⁡θ,mVMθ˙=Lcos⁡γV-Zsin⁡γV-mgcos⁡θ,-mVMcos⁡θψ˙V=Ysin⁡γV+Zcos⁡γV,x˙M=VMcos⁡θcos⁡ψV,y˙M=VMsin⁡θ,z˙M=-VMsin⁡θsin⁡ψV,where VM, θ, ΨV, and γV denote velocity, trajectory inclination angle, path angle, and bank angle of missile, respectively; xM, yM, and zM represent distance along x-, y-, z-axis; and m and g are mass of missile and acceleration of gravity, separately.Here the earth is assumed to be flat, and the model for the drag (D), the lift (L), and the side force (Z) can be formulated as(2)D=qSCD,L=qSCL,Z=qSCz,where q is dynamic pressure and S is reference area. Here the coefficients of aerodynamics CD, CL, and CZ are usually obtained by CFD or wind tunnel test.


### 2.2. Cost Function

The performance index used for minimization of missile miss-distance and control is(3)J=xtf2+∫t0tfαt2dt,where t denotes flight time and α is angle of attack (AOA). Here subscript 0 and f represent initial and terminal time, respectively.


### 2.3. Flight Path Constraints

The constraints corresponding to the flight path optimization for a given missile are as follows.The constraints on states and parameters control are(4)y∈ylower,yupper,θ∈θlower,θupper,where the subscript lower and upper represent lower and upper limit, respectively.These numbers are completely based on the flight envelope characteristics, representing the possible values for the maximum and minimum state variables.The terminal constraints of the missile are(5)Vtf≤Vf,θtf≥θf.Similarly, the nonlinear state inequality constraints are(6)αlower≤αt≤αupper,aylower≤ayt≤ayupper.For the purpose of the velocity control design, two subproblems are imperative to be solved. One is to obtain the optimal trajectory in the ideal condition. In consideration of this subproblem, the flight performance should be optimized to produce the trajectory and velocity profile satisfying the above requirements. The other one is to control the missile flying path along the optimal trajectory, in the same time to ensure the coherence of the missile speed and the ideal velocity profile; that is, the integration of the trajectory and velocity control should be considered.


## 2.1. Motion Equations

For the optimal flight velocity control problem under consideration, a mathematical model based on point mass dynamics can be used for determining the trajectory and velocity, yielding(1)mV˙M=-D-mgsin⁡θ,mVMθ˙=Lcos⁡γV-Zsin⁡γV-mgcos⁡θ,-mVMcos⁡θψ˙V=Ysin⁡γV+Zcos⁡γV,x˙M=VMcos⁡θcos⁡ψV,y˙M=VMsin⁡θ,z˙M=-VMsin⁡θsin⁡ψV,where VM, θ, ΨV, and γV denote velocity, trajectory inclination angle, path angle, and bank angle of missile, respectively; xM, yM, and zM represent distance along x-, y-, z-axis; and m and g are mass of missile and acceleration of gravity, separately.Here the earth is assumed to be flat, and the model for the drag (D), the lift (L), and the side force (Z) can be formulated as(2)D=qSCD,L=qSCL,Z=qSCz,where q is dynamic pressure and S is reference area. Here the coefficients of aerodynamics CD, CL, and CZ are usually obtained by CFD or wind tunnel test.


## 2.2. Cost Function

The performance index used for minimization of missile miss-distance and control is(3)J=xtf2+∫t0tfαt2dt,where t denotes flight time and α is angle of attack (AOA). Here subscript 0 and f represent initial and terminal time, respectively.


## 2.3. Flight Path Constraints

The constraints corresponding to the flight path optimization for a given missile are as follows.The constraints on states and parameters control are(4)y∈ylower,yupper,θ∈θlower,θupper,where the subscript lower and upper represent lower and upper limit, respectively.These numbers are completely based on the flight envelope characteristics, representing the possible values for the maximum and minimum state variables.The terminal constraints of the missile are(5)Vtf≤Vf,θtf≥θf.Similarly, the nonlinear state inequality constraints are(6)αlower≤αt≤αupper,aylower≤ayt≤ayupper.For the purpose of the velocity control design, two subproblems are imperative to be solved. One is to obtain the optimal trajectory in the ideal condition. In consideration of this subproblem, the flight performance should be optimized to produce the trajectory and velocity profile satisfying the above requirements. The other one is to control the missile flying path along the optimal trajectory, in the same time to ensure the coherence of the missile speed and the ideal velocity profile; that is, the integration of the trajectory and velocity control should be considered.


## 3. Velocity Control System Design

The process of the velocity control system design is conducted in Figure1. Primarily, the optimal trajectory and velocity profile are carried out as a criterion for meeting the index requirements. Then, in order to perform this velocity control, a virtual moving target moving along the critical trajectory is applied to guide the missile. In terms of this, the velocity control schematic and model can be built, and the velocity control law and control parameters will be implemented including the integration design for velocity control and trajectory control for a good performance of the coning motion missile system.Figure 1
The process of the velocity control system design.


### 3.1. Trajectory and Velocity Profile Optimization

With the advent of computers and evolution of modern theories of optimal control, the numerical computation techniques for optimal atmospheric trajectories have been an active research area since the early 1970s. The approaches on trajectory optimization are of two distinct categories such as direct method based on the mathematical programming and parameterization of state and control histories and the indirect method grounded on the solution of two-point boundary value problem (TPBVP) using optimal control principle [16–18]. Because the guess of the initial values of costate variables is random and there is a lack of physical implication, the TPBVP is difficult to solve especially when the optimal system is accompanied with high nonlinearity and multiple constraints (e.g., nonlinear trajectory optimization problem). The direct method has better convergence properties and thus can perform well with a poor initial guess [17]. From the mathematical point of view, the KKT conditions of NLP transformed by a discretization are equivalent to first-order necessary conditions of the original optimal control problem, which means the solution of the NLP is equivalent to that of the original optimal control problem.Therefore, the direct method is derived for the trajectory and the velocity profile optimization of the missile. Furthermore, the hp-adaptive Pseudospectral Method, as a kind of efficient direct method, is combining Legendre Pseudospectral Method [19, 20] and hp-adaptive method [21], which discretizes the state variables and control variables on a series of Legendre-Gauss-Lobatto (LGL) points. Through constructing a Lagrange interpolation polynomial to approach the state variables and the control variables by taking these discrete points as nodes, the optimal control problem is transformed into a nonlinear dynamic programming problem. Then, the constraints of differential equations can be altered into the form of algebra equations, and so do the integration item and the terminal status of the cost function; namely, they can be calculated by using Gauss-Lobatto integration and integrating the right function from initial status, respectively. When the calculation precision during some intervals fails to meet the requirement, the collocation numbers and order number of global interpolation polynomial should be adjusted adaptively in accordance with the hp-adaptive method.Generally, the discretization for the trajectory is conducted as follows. For transforming flight time [t0,tf] into the form of Legendre Pseudospectral Method [-1,1], the variable t is converted to(7)τ=2ttf-t0-tf+t0tf-t0.IfK is for the collocation numbers, the state vector and control vector will be approximately described by K Lagrange  polynomials Liτ considered as basis functions; that is,(8)xτ≈Xτ=∑i=1KLiτXτi,uτ≈Uτ=∑i=1KLiτUτi,where state vector x=[VMθψVxMyMzM] and control vector u=[α].Derive (7), and the approximate derivation of state vector is(9)x˙τk≈X˙τk=∑i=1KL˙iτkXτi.The constraint of dynamics equations thus is converted to the form of algebra constraint; that is,(10)∑i=1KDkiXτi-tf-t02fXτk=0,where Dki=L˙iτk means differentiation matrix of Legendre Pseudospectral Method.In terms of the above process of discretization, the interior-point method is applied to deal with the bounded constraints of the inequality such as (4), (5), and (6); that is, a barrier term μ is brought in to the objective function for completely avoiding this problem of bound constraints. Fortunately, C++ software pack (Ipopt) [22] of the interior-point method can be utilized to solve this discrete nonlinear programming (NLP). A design philosophy is described to enable the optimal generic design by using the hp-adaptive Pseudospectral Method, outlined briefly in Figure 2. Then, the schematic of algorithm to solve optimal control problem can be described as follows:(1)
Initialize a new mesh grid.(2)
Discretize the continuous optimization control problem through the Legendre Pseudospectral Method and transform it into NLP problem through computing LGR points (X(τi),U(τi)), weights Liτ, and differentiation matrices Dk.(3)
Solve the NLP problem by using the interior-point method.(4)
Update control variable and go back to Step (3), if the index J of cost function is not minimized.(5)
Quit the optimization process; otherwise update parameters and go back to Step (2), if the solution of state and path constraints meets error tolerance.Figure 2
Schematic of hp-adaptive pseudospectral algorithm to solve optimal control problem.


### 3.2. Trajectory Control by Tracking Virtual Moving Target

After obtaining the standard optimization trajectory and velocity control profile, the task shifts to pursuit of a design to make the missile fly along the ideal trajectory as well as the velocity profile. The conventional guidance method of the shaped trajectory can ensure the performance of the precision to the fixed target, which does not consider the requirement of the velocity control. Therefore, a virtual target proposed here is to afford a possibility of velocity control for the missile. Concretely, a virtual target is designed to move exactly like the ideal movement of the missile, that is, to move strictly along the standard optimized trajectoryAO⌢ with the same speed to the velocity profile. If a control design is developed to make the missile precisely chase the virtual target, the requirement of the velocity will be naturally satisfied for the missile flying along the ideal trajectory. Thus, in terms of the longitudinal plane of the kinematics described in Figure 3, all parameters of the position (xT,yT) and velocity VT including initial position (xT0,yT0), ideal initial velocity VT0, and anticipant terminal velocity VTf∗ can be obtained from the information of the optimized trajectory.Figure 3
Kinematical relationship of the missile and the virtual target.


### 3.3. Velocity Control Methodology

Figure4 depicts the total control flow of the coning motion control. It is known that the desired initial velocity VM0∗ (start from point A shown in Figure 3) can be conducted by using the above hp-adaptive pseudospectral trajectory optimization method for deciding whether the coning motion is need. Besides, in consideration of the influence of the control dynamics to the velocity, the coning motion will be made to descent of the flight velocity, when real velocity VM0 is more than 1.1 VM0∗.Figure 4
Control flow chart of the coning motion control.Obviously, the relative velocity and distance between the missile and the virtual target can be written asΔV=VM-VT and ΔR=RT-RM, respectively. Here, let ΔR be larger than zero when the target moves in front of the missile. If the missile motion satisfies the condition of ΔR>0 or VM<VT, the terminal velocity will surely meet the index constraint. Thus, the velocity control system can be built up in accordance with the relationship of the relative distance and speed between the missile and the virtual target (see Figure 5). Then, the velocity control law can be formulated as(11)axc=fΔV,ΔR=KRΔR-ΔRt+KVVT-VM,where ΔRt is the desired distance between the missile and the virtual target and KR and KV denote the velocity control gains.Figure 5
Velocity control structure of the missile.For the sake of convenience of analysis, it is good to suppose the anticipant distance between the missile and the virtual target to be0; that is, ΔRt=0. Therefore, the velocity control law will be simplified as(12)fΔV,ΔR=KRΔR+KVΔV.If the missile overlaps at the position of virtual target, for example,f(ΔV,ΔR)=0, ΔV=0, ΔR=0, the index of terminal velocity is satisfied naturally. Moreover, if f(ΔV,ΔR)<0, there are two circumstances: (1) the missile moves in back of the virtual target and (2) the missile is in front of the virtual target, but the speed magnitude is less than that of virtual target. Since the virtual target moves along the ideal flight trajectory all over, the circumstances of (1) and (2) will both satisfy the terminal velocity index. Additionally, the kinetic energy of the missile is unnecessarily large when f(ΔV,ΔR)>0, which will result in the missile deceleration by coning motion.Consequently, the cost functionf(ΔV,ΔR) determines whether the velocity control is required. For example, if f(ΔV,ΔR) is less than zero, the missile coning motion should be paused, and the proportional navigation guidance to the target should be switched on as well. Otherwise, if f(ΔV,ΔR) is greater than zero, the coning motion will be restarted.According to the control structure of the missile, dynamic inversion theory [23] can be used to design the control parameters KR and KV. The equation related to the velocity, shown in Figure 3, is expressed as(13)mV˙M=-qSCx-mgsin⁡θ.ReplacingV˙M with axc gives(14)axc=-qSCxm-gsin⁡θ.From the velocity control structure outlined in Figure4, transfer function VT to VM becomes(15)Ws=VMsVTs=KVs+KRs2+KVs+KR,where the 2nd-order characteristic equation is(16)s2+KVs+KR=0.Control gainsKV and KR are determined such that VM may become its desiring velocity VT. Apparently, from the linear system theory, the stability of the system is ascertained by the pole points, and the system is stable only when all pole points are located at the left of s plane.(1)
IfKV2≥4KR, since KV and KR are both over zero, the two characteristic roots of the system will be both minus real numbers, which can be written as(17)s1,2=-KV±KV2-4KR2.(2)
IfKV2<4KR, s1,2 changes to(18)s1,2=-KV±i4KR-KV22.Apparently, the nearer pole to the imaginary axis primarily determines the speed of output response of the above system. Therefore, if the poless1,2 have the same minus real part, that is, when KV2≤4KR, it is possible to obtain the fastest output response. Synchronously, the imaginary part of the characteristic roots determines the magnitude of the overshoot. In general, the relative damping value of system is optimally set to 2/2, and then the equation(19)4KR-KV2KV=22gives(20)KR=38KV2.Let VT(s)=1/s; substituting it into transfer function (15) yields(21)VMs=1s·KVs+KRs2+KVs+KR=1s-1s+0.5KV+0.5KVs+0.5KV2which can be described in the time-domain form by applying Laplace inverse transform; that is,(22)VMt=1+e-0.5KV·t0.5KV·t-1.Besides, by differentiating (22), it becomes(23)ddtVMt=12KVe-1/2KV·t2-12KVt.Let d/dtVM(t)=0, and the settling time of the system can be given as(24)Tp=4KV.In order to ensure the performance of the system, the settling timeTp is generally set to the triple of the time constant of the missile Td; namely,(25)Tp=3Td.By combining (20), (24), and (25), the solution to the velocity control law is(26)KV=43Td,KR=38KV2.When the missile does pure coning motion, the acceleration commands in longitudinal and lateral plane,aδyc and aδzc, are(27)aδyc=aδsin⁡ωt,aδzc=aδcos⁡ωt,where ω represents a circular frequency and aδ is the acceleration command corresponding to command of AOA δ in condition of the pure coning motion.In accordance with the velocity control law shown in (12), drag force command Xc can be acquired as(28)-Xc=maxc+mgsin⁡θwhich gives the lift force command(29)Yc=pXc=p·maxc+mgsin⁡θ,where p means lift-drag ratio. So(30)aδ=Ycm=pXcm=p·maxc+mgsin⁡θm≈pKVΔV+KRΔR.


### 3.4. Integrated Design for Velocity Control and Trajectory Control

On the other hand, in order to make the missile hit the target precisely, it is necessary to design an appropriate control law, under which the normal accelerationam can be formulated with the relative motion between the missile and the target as follows:(31)am=fr,ξ,where r denotes directive distance between missile and target and ξ is angle of line-of-sight (LOS).Traditional control laws are determined to approach the convergence of the distance and the LOS rate between the missile and the target. For instance, in terms of classic proportional navigation (PN) method, it is concluded that missile acceleration commandsam are in proportion to the rate of LOS ξ˙.From (27), the integrated guidance commands can be written as follows:(32)ayc=Δayc+aδsin⁡ωt,azc=Δazc+aδcos⁡ωt,where Δayc and Δazc denote the normal PN commands in longitudinal and lateral plane, respectively.As given in (32), if the velocity meets the anticipant requirement, that is, aδ=0, the missile will stop coning motion and will fly through the proportion guidance law so as to hit the target accurately. Otherwise, the missile will do the coning motion, which makes a continuous alternate changing between the angle of attack δsin⁡ωt and sideslip angle δcos⁡ωt, respectively. In an alternating cycle, the equivalent lift force produced by the alternate AOA is equal to zero, and the same to the side force. That is to say, only the drag force increases during the period of coning motion.Additionally, consistent with the normal control commandsΔayc and Δazc, the total AOA of the missile can be given as(33)Δ=α2+β2=Δα+δsin⁡ωt2+Δβ+δcos⁡ωt2=Δα2+Δβ2+hδ,Δα,Δβ+δ2,where β represents sideslip angle and h(δ,Δα,Δβ)=2δ(Δαsin⁡ωt+Δβcos⁡ωt).From (33), if the missile flies without coning motion, that is, δ=0, the total AOA will be simplified as(34)Δ=Δα2+Δβ2.Sinceh(δ,Δα,Δβ) changes periodically with time as shown in (33), it is possible that the missile with coning motion produces less total AOA than that without coning motion when the item h(δ,Δα,Δβ)+δ2<0. That is to say, coning motion control results in the missile’s drag decrease, especially when the weight of the missile gravity in the direction of the velocity is over the drag force; that is, -qSCx-mgsin⁡θ>0, and it is concluded that the speed of the missile will increase and not decrease. Consequently, to prevent this from happening, the inequation can be supposed as(35)Δαsin⁡ωt+Δβcos⁡ωt>0which gives(36)Δα2+Δβ2sin⁡ωt+ω0>0,where ω0 is equal to acos⁡(Δα/Δα2+Δβ2).In order to meet inequation (36) all along, circular frequency ω should satisfy(37)2πn≤ωt+ω0≤2πn+π,n∈E.So(38)2πn-ω0t≤ω≤2πn+π-ω0t.In view of the dynamic characteristics of the missile and the performance of the actuator, the bandwidth of the coining motion command must be less than that of the actuator. Otherwise, it is difficult for the missile guidance system to track the guidance command. Thus, the angular ratioω must be subject to ω<ωT; that is,(39)2πn+π-ω0t<ωTn=RoundωTt-π+ω02π,where ωT represents the bandwidth of the system.


## 3.1. Trajectory and Velocity Profile Optimization

With the advent of computers and evolution of modern theories of optimal control, the numerical computation techniques for optimal atmospheric trajectories have been an active research area since the early 1970s. The approaches on trajectory optimization are of two distinct categories such as direct method based on the mathematical programming and parameterization of state and control histories and the indirect method grounded on the solution of two-point boundary value problem (TPBVP) using optimal control principle [16–18]. Because the guess of the initial values of costate variables is random and there is a lack of physical implication, the TPBVP is difficult to solve especially when the optimal system is accompanied with high nonlinearity and multiple constraints (e.g., nonlinear trajectory optimization problem). The direct method has better convergence properties and thus can perform well with a poor initial guess [17]. From the mathematical point of view, the KKT conditions of NLP transformed by a discretization are equivalent to first-order necessary conditions of the original optimal control problem, which means the solution of the NLP is equivalent to that of the original optimal control problem.Therefore, the direct method is derived for the trajectory and the velocity profile optimization of the missile. Furthermore, the hp-adaptive Pseudospectral Method, as a kind of efficient direct method, is combining Legendre Pseudospectral Method [19, 20] and hp-adaptive method [21], which discretizes the state variables and control variables on a series of Legendre-Gauss-Lobatto (LGL) points. Through constructing a Lagrange interpolation polynomial to approach the state variables and the control variables by taking these discrete points as nodes, the optimal control problem is transformed into a nonlinear dynamic programming problem. Then, the constraints of differential equations can be altered into the form of algebra equations, and so do the integration item and the terminal status of the cost function; namely, they can be calculated by using Gauss-Lobatto integration and integrating the right function from initial status, respectively. When the calculation precision during some intervals fails to meet the requirement, the collocation numbers and order number of global interpolation polynomial should be adjusted adaptively in accordance with the hp-adaptive method.Generally, the discretization for the trajectory is conducted as follows. For transforming flight time [t0,tf] into the form of Legendre Pseudospectral Method [-1,1], the variable t is converted to(7)τ=2ttf-t0-tf+t0tf-t0.IfK is for the collocation numbers, the state vector and control vector will be approximately described by K Lagrange  polynomials Liτ considered as basis functions; that is,(8)xτ≈Xτ=∑i=1KLiτXτi,uτ≈Uτ=∑i=1KLiτUτi,where state vector x=[VMθψVxMyMzM] and control vector u=[α].Derive (7), and the approximate derivation of state vector is(9)x˙τk≈X˙τk=∑i=1KL˙iτkXτi.The constraint of dynamics equations thus is converted to the form of algebra constraint; that is,(10)∑i=1KDkiXτi-tf-t02fXτk=0,where Dki=L˙iτk means differentiation matrix of Legendre Pseudospectral Method.In terms of the above process of discretization, the interior-point method is applied to deal with the bounded constraints of the inequality such as (4), (5), and (6); that is, a barrier term μ is brought in to the objective function for completely avoiding this problem of bound constraints. Fortunately, C++ software pack (Ipopt) [22] of the interior-point method can be utilized to solve this discrete nonlinear programming (NLP). A design philosophy is described to enable the optimal generic design by using the hp-adaptive Pseudospectral Method, outlined briefly in Figure 2. Then, the schematic of algorithm to solve optimal control problem can be described as follows:(1)
Initialize a new mesh grid.(2)
Discretize the continuous optimization control problem through the Legendre Pseudospectral Method and transform it into NLP problem through computing LGR points (X(τi),U(τi)), weights Liτ, and differentiation matrices Dk.(3)
Solve the NLP problem by using the interior-point method.(4)
Update control variable and go back to Step (3), if the index J of cost function is not minimized.(5)
Quit the optimization process; otherwise update parameters and go back to Step (2), if the solution of state and path constraints meets error tolerance.Figure 2
Schematic of hp-adaptive pseudospectral algorithm to solve optimal control problem.


## 3.2. Trajectory Control by Tracking Virtual Moving Target

After obtaining the standard optimization trajectory and velocity control profile, the task shifts to pursuit of a design to make the missile fly along the ideal trajectory as well as the velocity profile. The conventional guidance method of the shaped trajectory can ensure the performance of the precision to the fixed target, which does not consider the requirement of the velocity control. Therefore, a virtual target proposed here is to afford a possibility of velocity control for the missile. Concretely, a virtual target is designed to move exactly like the ideal movement of the missile, that is, to move strictly along the standard optimized trajectoryAO⌢ with the same speed to the velocity profile. If a control design is developed to make the missile precisely chase the virtual target, the requirement of the velocity will be naturally satisfied for the missile flying along the ideal trajectory. Thus, in terms of the longitudinal plane of the kinematics described in Figure 3, all parameters of the position (xT,yT) and velocity VT including initial position (xT0,yT0), ideal initial velocity VT0, and anticipant terminal velocity VTf∗ can be obtained from the information of the optimized trajectory.Figure 3
Kinematical relationship of the missile and the virtual target.


## 3.3. Velocity Control Methodology

Figure4 depicts the total control flow of the coning motion control. It is known that the desired initial velocity VM0∗ (start from point A shown in Figure 3) can be conducted by using the above hp-adaptive pseudospectral trajectory optimization method for deciding whether the coning motion is need. Besides, in consideration of the influence of the control dynamics to the velocity, the coning motion will be made to descent of the flight velocity, when real velocity VM0 is more than 1.1 VM0∗.Figure 4
Control flow chart of the coning motion control.Obviously, the relative velocity and distance between the missile and the virtual target can be written asΔV=VM-VT and ΔR=RT-RM, respectively. Here, let ΔR be larger than zero when the target moves in front of the missile. If the missile motion satisfies the condition of ΔR>0 or VM<VT, the terminal velocity will surely meet the index constraint. Thus, the velocity control system can be built up in accordance with the relationship of the relative distance and speed between the missile and the virtual target (see Figure 5). Then, the velocity control law can be formulated as(11)axc=fΔV,ΔR=KRΔR-ΔRt+KVVT-VM,where ΔRt is the desired distance between the missile and the virtual target and KR and KV denote the velocity control gains.Figure 5
Velocity control structure of the missile.For the sake of convenience of analysis, it is good to suppose the anticipant distance between the missile and the virtual target to be0; that is, ΔRt=0. Therefore, the velocity control law will be simplified as(12)fΔV,ΔR=KRΔR+KVΔV.If the missile overlaps at the position of virtual target, for example,f(ΔV,ΔR)=0, ΔV=0, ΔR=0, the index of terminal velocity is satisfied naturally. Moreover, if f(ΔV,ΔR)<0, there are two circumstances: (1) the missile moves in back of the virtual target and (2) the missile is in front of the virtual target, but the speed magnitude is less than that of virtual target. Since the virtual target moves along the ideal flight trajectory all over, the circumstances of (1) and (2) will both satisfy the terminal velocity index. Additionally, the kinetic energy of the missile is unnecessarily large when f(ΔV,ΔR)>0, which will result in the missile deceleration by coning motion.Consequently, the cost functionf(ΔV,ΔR) determines whether the velocity control is required. For example, if f(ΔV,ΔR) is less than zero, the missile coning motion should be paused, and the proportional navigation guidance to the target should be switched on as well. Otherwise, if f(ΔV,ΔR) is greater than zero, the coning motion will be restarted.According to the control structure of the missile, dynamic inversion theory [23] can be used to design the control parameters KR and KV. The equation related to the velocity, shown in Figure 3, is expressed as(13)mV˙M=-qSCx-mgsin⁡θ.ReplacingV˙M with axc gives(14)axc=-qSCxm-gsin⁡θ.From the velocity control structure outlined in Figure4, transfer function VT to VM becomes(15)Ws=VMsVTs=KVs+KRs2+KVs+KR,where the 2nd-order characteristic equation is(16)s2+KVs+KR=0.Control gainsKV and KR are determined such that VM may become its desiring velocity VT. Apparently, from the linear system theory, the stability of the system is ascertained by the pole points, and the system is stable only when all pole points are located at the left of s plane.(1)
IfKV2≥4KR, since KV and KR are both over zero, the two characteristic roots of the system will be both minus real numbers, which can be written as(17)s1,2=-KV±KV2-4KR2.(2)
IfKV2<4KR, s1,2 changes to(18)s1,2=-KV±i4KR-KV22.Apparently, the nearer pole to the imaginary axis primarily determines the speed of output response of the above system. Therefore, if the poless1,2 have the same minus real part, that is, when KV2≤4KR, it is possible to obtain the fastest output response. Synchronously, the imaginary part of the characteristic roots determines the magnitude of the overshoot. In general, the relative damping value of system is optimally set to 2/2, and then the equation(19)4KR-KV2KV=22gives(20)KR=38KV2.Let VT(s)=1/s; substituting it into transfer function (15) yields(21)VMs=1s·KVs+KRs2+KVs+KR=1s-1s+0.5KV+0.5KVs+0.5KV2which can be described in the time-domain form by applying Laplace inverse transform; that is,(22)VMt=1+e-0.5KV·t0.5KV·t-1.Besides, by differentiating (22), it becomes(23)ddtVMt=12KVe-1/2KV·t2-12KVt.Let d/dtVM(t)=0, and the settling time of the system can be given as(24)Tp=4KV.In order to ensure the performance of the system, the settling timeTp is generally set to the triple of the time constant of the missile Td; namely,(25)Tp=3Td.By combining (20), (24), and (25), the solution to the velocity control law is(26)KV=43Td,KR=38KV2.When the missile does pure coning motion, the acceleration commands in longitudinal and lateral plane,aδyc and aδzc, are(27)aδyc=aδsin⁡ωt,aδzc=aδcos⁡ωt,where ω represents a circular frequency and aδ is the acceleration command corresponding to command of AOA δ in condition of the pure coning motion.In accordance with the velocity control law shown in (12), drag force command Xc can be acquired as(28)-Xc=maxc+mgsin⁡θwhich gives the lift force command(29)Yc=pXc=p·maxc+mgsin⁡θ,where p means lift-drag ratio. So(30)aδ=Ycm=pXcm=p·maxc+mgsin⁡θm≈pKVΔV+KRΔR.


## 3.4. Integrated Design for Velocity Control and Trajectory Control

On the other hand, in order to make the missile hit the target precisely, it is necessary to design an appropriate control law, under which the normal accelerationam can be formulated with the relative motion between the missile and the target as follows:(31)am=fr,ξ,where r denotes directive distance between missile and target and ξ is angle of line-of-sight (LOS).Traditional control laws are determined to approach the convergence of the distance and the LOS rate between the missile and the target. For instance, in terms of classic proportional navigation (PN) method, it is concluded that missile acceleration commandsam are in proportion to the rate of LOS ξ˙.From (27), the integrated guidance commands can be written as follows:(32)ayc=Δayc+aδsin⁡ωt,azc=Δazc+aδcos⁡ωt,where Δayc and Δazc denote the normal PN commands in longitudinal and lateral plane, respectively.As given in (32), if the velocity meets the anticipant requirement, that is, aδ=0, the missile will stop coning motion and will fly through the proportion guidance law so as to hit the target accurately. Otherwise, the missile will do the coning motion, which makes a continuous alternate changing between the angle of attack δsin⁡ωt and sideslip angle δcos⁡ωt, respectively. In an alternating cycle, the equivalent lift force produced by the alternate AOA is equal to zero, and the same to the side force. That is to say, only the drag force increases during the period of coning motion.Additionally, consistent with the normal control commandsΔayc and Δazc, the total AOA of the missile can be given as(33)Δ=α2+β2=Δα+δsin⁡ωt2+Δβ+δcos⁡ωt2=Δα2+Δβ2+hδ,Δα,Δβ+δ2,where β represents sideslip angle and h(δ,Δα,Δβ)=2δ(Δαsin⁡ωt+Δβcos⁡ωt).From (33), if the missile flies without coning motion, that is, δ=0, the total AOA will be simplified as(34)Δ=Δα2+Δβ2.Sinceh(δ,Δα,Δβ) changes periodically with time as shown in (33), it is possible that the missile with coning motion produces less total AOA than that without coning motion when the item h(δ,Δα,Δβ)+δ2<0. That is to say, coning motion control results in the missile’s drag decrease, especially when the weight of the missile gravity in the direction of the velocity is over the drag force; that is, -qSCx-mgsin⁡θ>0, and it is concluded that the speed of the missile will increase and not decrease. Consequently, to prevent this from happening, the inequation can be supposed as(35)Δαsin⁡ωt+Δβcos⁡ωt>0which gives(36)Δα2+Δβ2sin⁡ωt+ω0>0,where ω0 is equal to acos⁡(Δα/Δα2+Δβ2).In order to meet inequation (36) all along, circular frequency ω should satisfy(37)2πn≤ωt+ω0≤2πn+π,n∈E.So(38)2πn-ω0t≤ω≤2πn+π-ω0t.In view of the dynamic characteristics of the missile and the performance of the actuator, the bandwidth of the coining motion command must be less than that of the actuator. Otherwise, it is difficult for the missile guidance system to track the guidance command. Thus, the angular ratioω must be subject to ω<ωT; that is,(39)2πn+π-ω0t<ωTn=RoundωTt-π+ω02π,where ωT represents the bandwidth of the system.


## 4. Simulation Example

In this section, a simulation example of the velocity control system by using coning motion technique is demonstrated for a missile, whose indexes include terminal speed, terminal path angle, and miss distance. In order to exhibit the velocity control performance, the conventional PN mode without the coning motion and particle swarm optimization (PSO) approach are involved for comparing study. The initial conditions that the missile makes coning motion are set as follows: the initial values ofy, θ, V, and ΨV are set to be x(0)=-19km, y(0)=10km, θ(0)=-20deg, V(0)=500m/s, and ΨV(0)=0deg. The initial position of the virtual target is selected as (-7292,9329,0)m besides the total flight time tT 63.8 s. The flight constraints are set as θlow=-30°, θtop=30°, ylow=0km, ytop=10km, Vf=250m/s, θf=-25°, αdown=-20deg, αup=20deg, aydown=-50m/s2,   and ayup=50m/s2. Additionally, the design process of PSO trajectory is referred to [24], and the parameters of the adopted PSO method are m=50, n=2, and c1=c2=1.8. Inertial weight w is set to be reduced linearly, from 1.0 to 0.4, and the search will be terminated if the fitness value is less than 0.4 or the number of iteration reaches 150.According to the above integrated design method for velocity control and trajectory control, the simulation results are illustrated in Figures6–13. Comparing the results of the optimal method proposed in this paper with those of the PSO approach, it is concluded that the optimization solutions of two approaches are almost the same, which implies the presented direct discretization method has sufficient precision of optimization stratifying the requirements of the constraints as exhibited in Table 1.Table 1
Simulation results produced by different guidance laws.
Guidance with velocity control
Guidance without velocity control
Optimization trajectory profile
PSO trajectory profile
Final velocity (m/s)
254.09
281.36
255.0
255.6
Final flying path angle (deg)
−15.65
−20.18
−21.46
−20.93
Miss distance (m)
2.2
1.6
0.6
0.8Figure 6
3D curve of the flight trajectories with different modes.Figure 7
History curve of the coning motion.Figure 8
Variation of the nutation angle.Figure 9
Response of the velocity in different condition.Figure 10
Variation of the angle of attack.Figure 11
Variation of the sideslip angle.Figure 12
Curve of the trajectory inclination angle.Figure 13
Curve of the path angle versus range.Additionally, as the result shown in Figure6, the velocity control system has got the same normal control performance to the conventional proportional navigation scheme. From Figures 7-8 where points A and B denote the start and end of the coning motion, respectively, the control system improves the magnitude amount of induced AOA, which makes the missile’s drag increase. Therefore, the speed of the missile decreases apparently in this period, as illustrated in Figure 9. Meanwhile, Figures 10-11 indicate that the coning motion control suspends at the range of 13 km and switches into the pure proportion guidance mode, which suggests that the missile velocity is gradually moving to the index of requirement because the magnitudes of the attack and sideslip angles are both getting smaller with the missile approaching to the virtual target. According to the path angle curves outlined in Figures 12-13, it is helpful to suspend the coning motion control before hitting the target because the coning motion reduces the tracking precision to the ideal trajectory. Finally, as shown in Table 1, both guidance modes meet the requirements of high guidance precision and terminal path angle, but the coning motion control system produces lower amount of terminal velocity, which is in line with the desired velocity index.


## 5. Conclusions

In this paper, a novel coning motion based guidance method is proposed for the velocity control of missile system with terminal velocity constraint. An ideal velocity profile is presented by using nonlinear programming method, in which a virtual moving target is put forward for missile to chase so as to perform the velocity control. Moreover, the dynamic inversion theory is applied to design the velocity control system parameters. To show the effectiveness of the control law, comparative simulations of a missile, in which pure proportional navigation mode has been taken into consideration, are provided. Simulation results demonstrate that the proposed velocity control law has a good performance providing reference to the design for velocity control system of the missile. The extensions of the presented coning motion control laws to a spinning missile and the inclusion of the dynamics in the velocity control system are possible areas of further research.



---

*Source: 716547-2015-12-06.xml*