A novel quasi-dynamic guidance law for a dynamic dual-spin projectile with non-conventional,asymmetric roll constraints

Abstract

A novel quasi-dynamic guidance law (QDGL) is presented for a dual-spin projectile (DSP) with unconventional constraints on roll direction. A 7 degree-of-freedom (DOF) dynamic model is established and the projectile operational mechanism is presented with a description of how it is used to enact control. The QDGL is presented and a parametric study is conducted to show how the QDGL parameters affect the system response. A procedure of using batches of Monte Carlo simulations is described, to numerically compare the system response with different QDGL configurations. A genetic algorithm is then used to optimise both the innate system parameters and PID controller gains. The disturbance rejection capabilities of the optimal QDGL are then evaluated along with the performance against different target profiles. It was found that the GA optimised QDGL is able to provide satisfactory control capabilities against static and dynamic targets.

Keywords

Guidance law dual-spin projectile genetic algorithm

Introduction

The calibre of guided weapons is getting ever smaller to meet the evolving needs of modern engagement scenarios. This is in turn driving a reduction of critical subsystem volumes of guidance and control hardware. It is thus pertinent to develop a novel method of control so as to either minimise the necessary subsystems occupying a volume or reduce the demand on subsystems/materials which are included in the projectile. This would allow a reduction in calibre with currently available technology without having to advance the capabilities of any specific subsystem/material.

Dual-spin projectiles (DSPs), such as STARSTREAK,¹ are becoming more prevalent in today’s military arsenals as they enable a wide scope of target engagement profiles. The dual-spin configuration slows the roll rate of the forward section to a point where the response rate of the actuators is sufficiently high compared to their roll rate such that effective control can be enacted. The aft section stays at the high roll rate and the projectile thus maintains gyroscopic stability. Guidance modules in the form of course corrected fuses are also being retrofitted onto conventional munitions, such as the Orbital ATK Armament Systems’ M1156 Precision Guidance Kit² and BAE’s Silver Bullet.³ In these larger weapon systems (155 mm), the roll rate of the projectiles is relatively low compared to smaller calibres. For smaller projectile calibres with high roll rates, it may not prove feasible to mitigate the high roll rates to apply conventional control methods.

Some prevalent examples of conventional guidance laws (GLs) used in projectiles include proportional navigation (PN), proportional derivative (PD) and sliding mode control (SMC).⁴ Conventional SMC or variations thereof have already considered constraints such as autopilot lag and actuator fault,⁵ accelerator saturation,⁶ and modelling uncertainty in missile/target dynamics. Impact angle is often the considered an aspect for control in a GL. In addition to purely controlling the impact angle,⁷ secondary constraints have also been placed on trajectory time,⁸ field-of-view,^9–11 and manoeuvrability.¹² External uncertainties and missile jerk have also been considered.¹³

A modified SMC GL was used to improve the chattering, miss-distance and finite time over conventional SMC and PN methods.¹⁴ The validity of the PN-like LOS GL has been investigated for a three body (two aircraft, one missile) system where the launch platform is also moving.¹⁵ A novel variation, called ‘airborne-CLOS’ utilises two separate LOS rates with one gain to control the three body problem.¹⁶

A novel GL has been created utilising virtual targets for impact angle and burst height constraints.¹⁷ A polar GL has been investigated which controls a missile based on the polar radius and angle of the target from the missile.¹⁸ An expanded 2D PD GL was created for a skid-to-turn command to LOS anti-tank guided missile, which builds upon classic PD, with the objective of eliminating a spiral trajectory which is an artefact of PD GLs.¹⁹ A proposed method uses a weighted zero-effort-miss (ZEM) to shape the actual ZEM, presenting as a PN GL with an extra time varying gain.

A few publications specifically pertain to the guidance of DSPs. Iterative impact point prediction has been used to create a GL for a DSP with control force imparted by fixed canards²⁰. A modified form of projectile linear theory is used to predict where the projectile will land and make the necessary corrections to the control system. Proportional navigation has been used in the GL of a dual-spin mortar during the ascent and descent phase.²¹ The results of the GL were validated with hardware-in-the-loop testing and Monte Carlo simulations.

It is a common practice to neglect gravitational forces when creating a GL.^22–25 Once an idealistic GL has been created, constraints can then be placed on the model which reflect conditions present in chosen real-world systems. GLs are often described from the perspective of the YZ plane (shown in Figure 1) also known as the ‘picture plane’.^7,18,19 Some of the cited literature use kinematic models for the derivation and validation of the control law,^{8,12,14,15,17} while this paper uses a dynamic projectile model. It is common to test the GL using arbitrary model parameters to facilitate more efficient and reliable interpretation of the results.^18,19

Figure 1.

Dual-spin projectile with fixed control force F_c.

Conventional projectile control utilises control surfaces which are able to adjust the roll angle of a projectile as well as the magnitude of the control force. The conventional guidance strategy is to roll the projectile to align the controllable pitch axis with the desired direction, and then increase the force by actuating the control surfaces which results in lateral movement. Dual-spin projectiles use a similar method with the addition of a coaxial motor to assist the correction of the forward section roll angle.²⁶ This paper proposes a quasi-dynamic GL (QDGL) for a DSP, with a fixed roll direction as well and a fixed magnitude control force; however, the phase of the force can change, hence the quasi nature. Control is enacted by adjusting the roll rate of the control force, slowing it down through certain roll angles to bias the force in the desired direction. There currently exists no literature which describes a GL for a DSP with roll-direction, roll-rate and control-force magnitude constraints.

The section Projectile Dynamics describes the 7 degree-of-freedom (DOF) dynamic model of the projectile used in simulations. The section Quasi-dynamic guidance law formulation introduces the projectile design and describes how control is enacted with the asymmetric roll constraints. The QDGL and associated parameters are introduced, with a brief parametric investigation showing their effect on the system response. The section Monte Carlo Procedure and Normalised Errors presents a Monte Carlo procedure which is used to numerically compare the effect that different QDGL configurations have on the system response. A genetic algorithm is then run to optimise the QDGL parameters as well as the gains of a PID controller. The system response of the optimised QDGL is evaluated against various disturbances and target profiles. The final section gives a summary of the paper and the key findings.

Projectile dynamics

This section describes the seven DOF dynamic model of the DSP, which has previously been used to investigate DSPs^21,27–29 and is a derivative of the well-established model for a conventional projectile used by McCoy.³⁰ Subscript F denotes the forward section and subscript A denotes the aft section. The assumption is made that the total centre of mass (COM) coincides with the aft COM, that is, the mass of the forward part is small with respect to that of the aft part and the nose moment of inertia Ixx,F is small compared to the aft one Ixx,A. The forces and moments the projectile is subject to are represented by an aerodynamic coefficient. The whole body longitudinal C_A0, transverse C_Y0 and normal C_N0 coefficients represent the combined effect of these individual forces, and are shown in Figure 2 for non-zero angles of attack α.

Figure 2.

Projectile coefficients when α ≠ 0.

The non-linear kinematic translational and rotational equations are given by equations (1) and (2), respectively

[\begin{matrix} {\dot{x}}_{e} \\ {\dot{y}}_{e} \\ {\dot{z}}_{e} \end{matrix}] = [\begin{matrix} \cos (θ) \cos (ψ) & - \sin (ψ) & \sin (θ) \cos (ψ) \\ \cos (θ) \sin (ψ) & \cos (ψ) & \sin (θ) \sin (ψ) \\ - \sin (θ) & 0 & \cos (θ) \end{matrix}] [\begin{matrix} u \\ v \\ w \end{matrix}]

(1)

[\begin{matrix} {\dot{ϕ}}_{F} \\ {\dot{ϕ}}_{A} \\ \dot{θ} \\ \dot{ψ} \end{matrix}] = [\begin{matrix} 1 & 0 & 0 & \tan (θ) \\ 0 & 1 & 0 & \tan (θ) \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & \sec (θ) \end{matrix}] [\begin{matrix} p_{F} \\ p_{A} \\ q \\ r \end{matrix}]

(2)

Accordingly, the dynamic translational is shown in equation (3)

[\begin{matrix} \dot{u} \\ \dot{v} \\ \dot{w} \end{matrix}] = \frac{1}{m} \begin{matrix} F_{X} \\ F_{Y} \\ F_{Z} \end{matrix}] - [\begin{matrix} 0 & - r & q \\ r & 0 & r \tan (θ) \\ - q & - r \tan (θ) & 0 \end{matrix}] \cdot [\begin{matrix} u \\ v \\ w \end{matrix}]

(3)

while the dynamic rotational equation is shown in equation (4)

[\begin{matrix} {\dot{p}}_{F} \\ {\dot{p}}_{A} \\ \dot{q} \\ \dot{r} \end{matrix}] = [\begin{matrix} 0 \\ 0 \\ - (\frac{I_{x x, A}}{I_{y y}}) p_{A} r - r^{2} \tan (θ) \\ (\frac{I_{x x, A}}{I_{y y}}) p_{A} q + q r \tan (θ) \end{matrix}] + [\begin{matrix} I_{x x, F}^{- 1} & 0 & 0 & 0 \\ 0 & I_{x x, A}^{- 1} & 0 & 0 \\ 0 & 0 & I_{y y}^{- 1} & 0 \\ 0 & 0 & 0 & I_{y y}^{- 1} \end{matrix}] \cdot [\begin{matrix} m_{1, F} \\ m_{1, A} \\ m_{m} \\ m_{n} \end{matrix}]

(4)

Here, the forces ${[F_{X} F_{Y} F_{Z}]}^{T}$ and moments ${[m_{1, F} m_{1, A} m_{m} m_{n}]}^{T}$ are given by equations (5) and (6), respectively

[\begin{matrix} F_{X} \\ F_{Y} \\ F_{Z} \end{matrix}] = \bar{q} S {[\begin{matrix} - C_{A 0} \\ C_{Y 0} \\ - C_{N 0} \end{matrix}] + \frac{p_{A} d}{V} [\begin{matrix} 0 \\ C_{Y p} \\ - C_{N p} \end{matrix}] + [\begin{matrix} 0 \\ C_{Y δ} δ_{y} \\ - C_{N δ} δ_{z} \end{matrix}]} + m g [\begin{matrix} - \sin (θ) \\ 0 \\ \cos (θ) \end{matrix}]

(5)

[\begin{matrix} m_{1, F} \\ m_{1, A} \\ m_{m} \\ m_{n} \end{matrix}] = \bar{q} S d {[\begin{matrix} 0 \\ 0 \\ C_{m 0} \\ - C_{n 0} \end{matrix}] + \frac{d}{V} [\begin{matrix} 0 \\ C_{l p} \\ C_{m q} \\ C_{n r} \end{matrix}] + \frac{p_{A} d}{V} [\begin{matrix} 0 \\ 0 \\ C_{m p} \\ C_{n p} \end{matrix}] + [\begin{matrix} 0 \\ 0 \\ C_{m δ} δ_{z} \\ C_{n δ} δ_{y} \end{matrix}]} + [\begin{matrix} τ_{M} + τ_{F, A} \\ - τ_{F, A} \\ 0 \\ 0 \end{matrix}]

(6)

Here, $ℳ$ is the Mach number, τ_M is the motor torque which controls the roll angle of the forward guidance section (ϕ_F) and τ_FA is the frictional torque which counteracts this. The frictional torque can be modelled as a hydrodynamic and roller bearing in combination where k_s and k_v are the static and viscous friction coefficients, respectively

τ_{F, A} = \frac{1}{2} ρ S d C_{A 0} (ℳ, α, β) sgn (p_{A} - p_{F}) (k_{s} + k_{v} | p_{A} - p_{F} |)

(7)

Quasi-dynamic guidance law formulation

This section describes the control method of the DSP design with unconventional roll constraints. The QDGL is formulated and the resulting system response is shown. A brief parametric study is conducted to illustrate how the QDGL parameters affect the system response.

Most guided projectiles have fins which can roll the projectile and induce a variable control force along a pitching axis. If lateral deflection is required, the projectile adjusts its roll angle such that the axis of the control force is parallel to the direction of required travel, whereby it increases the control force and therefore imparts an acceleration. Figure 3 shows the DSP design; a control force F_c is produced by aerodynamic lifting surfaces on the forward section. At launch, the aft section engages with the rifling which accelerates the roll rate, while the forward section producing F_c remains de-spun. During flight, the two sections will reach an equilibrium through the bearing torque τF,A, and the forward section will have a relatively slow roll rate, ω₀. If τF,A is increased during flight (e.g. by means of a brake) then there will be new equilibrium where the forward section has a higher roll rate, ω₁.

Figure 3.

Earth axis perspective of the picture plane and control force F_c rotating at rate ω₁ through angle ϕ

Figure 3 shows the key parameters of the control method, as well as the YZ plane, referred to as the picture plane. The constant magnitude F_c moves through a roll angle ϕ with rate ω₀ or ω₁, where ω₀ < ω₁. The ϕ ∈ [0, 2π] describes the roll orientation of F_c with respect to the normal axis. It sweeps in the negative mathematical direction, since most conventional projectiles have a right hand twist. The novel guidance strategy proposed herein uses a fixed magnitude F_c rolling at speed ω₁. The roll rate is slowed to ω₀ through favourable roll angles when F_c is aligned with the desired correction axis, and then accelerated back to ω₁ through the remaining unfavourable roll angles. The act of slowing F_c when sweeping through favourable roll angles is henceforth referred to as a ‘bias’. All measurements and symbols henceforth are given in the ‘picture plane’ reference frame unless explicitly stated.

The integral of Newton’s second law relates the impulse of an object, J, to its change in velocity Δv

J |_{Δ t} = m Δ v |_{Δ t}

(8)

Here, the mass is an assumed constant since there are no on-board resources being consumed. A generalised decomposition of F_c onto any orthonormal axis i, j in the YZ plane has the corresponding forces F_i, F_j. Let the desired decomposition axis i be an angle ϕ_B from the normal axis $\hat{z}$ (where ϕ = 0). Let ϕ_i be a particular angle between F_c and the arbitrary decomposition axis i. Let ϕ_a be the angle through which F_c sweeps at a given rate ω such that the sweep begins at the angle (ϕ_B − ϕ_a) and ends at ϕ_B.

The range of angles during which F_c is slowed is defined as the bias angle. Let the midpoint of the bias angle coincide with decomposition axis i, such that the symmetrical angle on either side of the midpoint is ϕ_a. The bias angle thus starts at (ϕ_B − ϕ_a) and ends at (ϕ_B + ϕ_a) with a midpoint of ϕ_B. This is shown in Figure 4. However, F_c will continue to rotate through the rest of the angle ϕ eventually sweeping another angular range (ϕ_B + π) ± ϕ_a (wrapped so ϕ ∈ [0, 2π]). During this time, the resulting change in velocity is directed along the negative i^th axis.

Figure 4.

Bias manoeuvre of size ϕ_a centred about ϕ_B.

We define ΔV as the total change in velocity of one whole roll rotation in sweeping through equal but opposing angles of size 2ϕ_a, at different rates ω₀ and ω₁. Assuming F_c, m and ω as constants, it can be shown from equation (8) that

Δ V = \frac{2 F_{c}}{m} \sin (ϕ_{a}) (\frac{ω_{0} - ω_{1}}{ω_{0} ω_{1}})

(9)

The maximum bias angle is half of a roll rotation, ϕ_a,max = π/2. The maximum ΔV per rotation is thus given by

Δ V_{max} = {Δ V |}_{ϕ_{a} = π / 2}

(10)

which is evaluated for a given system. Table 1 shows a comparison between the idealised parameters which are used during analysis of the QDGL, and values that can be expected for real systems. Lloyd and Brown³¹ investigated the maximum lateral control force F_c which can be applied to the nose of a projectile before flight instability occurs. For a 15 kg, 105 mm projectile, not more than 40N (0.34 g) could be applied. Li et al.²⁹ continued this work and concluded that the same projectile would remain stable provided F_c ∈ [ − 35.48, 58.33]. The provided F_c is within the limits defined by these frameworks; projectile stability is ensured by satisfying the gyroscopic stability criterion.³² When designing such a projectile, the gyroscopic and dynamic stability framework for DSPs³³ or other methods³⁴ can be used.

Table 1.

Comparison of simulation parameters used against real-world projectiles.

Parameter	Idealised system	105 mm	155 mm projectile
ω₀ [rad.s⁻¹]	π/2	N/A	N/A
ω₁ [rad.s⁻¹]	2π	1050²⁹	1668³⁵
m [kg]	1	15²⁹	42.798³⁶
F_c [N]	1	[ − 35.48, 58.33]²⁹	20³⁵

Since the purpose of this paper is to test whether the uni-rotational, fixed magnitude control force can be used to guide a projectile in the 2D plane, all forces and moments except the control force are neglected. Additionally, there is no transient between the fast and slow oscillations, and the switching is instantaneous.

By design, the QDGL calculates a desired change in speed when ϕ = 0, and then calculates the bias angles from equation (9). The projectile will then continue to roll, whereby the actuator will slow the roll down if the current roll angle lies within the bias range previously calculated. In practice, the desired speed change and resulting bias angles are calculated when ϕ lies in a small range, ϕ ∈ [0, 0.001], to account for the machine computation inaccuracy. While this calculation could be conducted and updated continuously, the relative speeds would have to be transformed to the ϕ = 0 reference frame which adds another layer of computational complexity. In addition, this discrete computation of speeds at the beginning of each rotation accommodates the bandwidth of hardware with respect to the roll rate of the projectile.

The current relative velocity of the projectile to target is the difference between the projectile and target velocity, V_R = V − V_T, or in full

V_{R} = [\begin{matrix} u_{R} \\ v_{R} \end{matrix}] = [\begin{matrix} u - u_{T} \\ v - v_{T} \end{matrix}]

(11)

N.B. The projectile having [u₀, v₀] = [0, 0] and undergoing consecutive unbiased rotations does not result in a circular trajectory; instead a semi-circular trajectory would result. To achieve a circular trajectory in the resting state, the horizontal velocity at the beginning of the bias calculation must assume the control force has already rotated through one quarter rotation. Taking this into consideration, we define V_DR0 as the ΔV correction necessary to bring the projectile to a stable circular orbit relative to the target, including the current relative velocity

V_{D R 0} = [u_{R} + Δ V |_{ϕ = π / 4} v_{R}]

(12)

This only allows the autopilot to bring the projectile to relative rest, the desired closing speed V_PT(d) describes the chosen approach speed as a function of d. The total demanded velocity change from the velocity autopilot V_Dem is then a linear combination of the necessary relative speed correction to bring the system to an orbit, V_DR0, and the closing velocity V_PT(d) dictated by the QDGL is

V_{Dem} = V_{D R 0} + V_{P T} (d)

(13)

V_PT(d) must only demand speeds which can be delivered by the actuation mechanism, given that ΔV can never exceed ΔV_max. Let the function V_lim(d) be the maximum relative speed the projectile can have at a distance d ≥ 0, such that it is still able to decelerate in time to be at relative rest when d = 0. This function can be calculated by starting with a stationary projectile and applying consecutive ΔV_max biases, since the process is reversible. From the rates given by the idealised system parameters in Table 1, a ΔV_max = 0.954 9 ms⁻¹ bias is enacted by the projectile in 2.5s. An effective acceleration value, $\tilde{a}$ , is measured from simulations for consecutive ΔV_max biases. Using this, it can be shown that

V_{Lim} (d) = {(2 \tilde{a} d)}^{1 / 2}

(14)

Since the function V_PT(d) is calculated when ϕ = 0 at a particular distance d₁, the desired ΔV will not be achieved until after the bias manoeuvre has been executed one full rotation later. Hence, the process is discontinuous. By this point, the projectile will have moved to some new distance d₂ under its residual velocity. Figure 5 shows the delay of the system response when the target speed is set to be the limit of the system, for example,. V_PT(d) = V_lim(d). The data points indicate specific positions where ϕ = 0, triggering the calculation of the bias angles. Figure 5 illustrates how the value of V_lim(d) demanded at a specific point is achieved at the next calculation point but after the delay caused by the roll rotation. This delay causes the system to exceed V_lim(d), resulting in an overshoot. To account for the delay, the demanded speed is modified by a factor ξ which ensures the relative speed never exceeds V_lim(d). The delay does not directly scale with distance but rather with V_PT(d) as it is the result of dynamic system evolution. Hence, the closing speed function is written as

V_{P T} (d) = V_{lim} (d) - ξ, ξ \in ℝ \geq 0

(15)

where ξ is a constant to be optimised. The result of this modification is also shown in Figure 5, where the adapted system response never exceeds the limit, for example,. V_PT(d) = V_lim(d) − ξ | V_R(d) > V_lim(d).

Figure 5.

Effect of ξ to prevent system response exceeding V_lim(d).

By including ξ, V_PT(d) is not properly defined when d ≤ d₁ where d₁| (V_PT(d₁) − ξ = 0). From equation (14), this boundary is

d_{1} = \frac{ξ^{2}}{2 \tilde{a}}

(16)

As such, the function V_lim(d) − ξ is only valid for d ∈ (d₁, ∞] and must be defined by other means for d ∈ [0, d₁]. Firstly, the distance d₂ is chosen to represent the desired level of precision for the projectile. The projectile will remain relatively stationary within this threshold, when it is on course to hit the target, so V_PT(d) = 0 ∀ d ∈ [0, d₂] where $d_{2} \in ℝ \geq 0$ . Secondly, a linear region is defined where the projectile moves at a low constant velocity toward the target in anticipation of either stopping or accelerating by a larger amount. This is represented by V_PT(d) = V_k ∀ d ∈ [d₂, d₁] where d₂ < d₁ and $V_{k} \in ℝ \geq 0$ is a constant speed. In its entirety, the function V_PT(d), and thus the QDGL, for d ∈ [0, ∞] is given by

V_{P T} (d) = [\begin{matrix} V_{\lim} (d) - ξ \\ V_{k} \\ 0 \end{matrix}] \forall d \in {\begin{array}{l} d_{1} \leq d \\ d_{2} \leq d < d_{1} \\ 0 \leq d < d_{2} \end{array}

(17)

Monte Carlo Procedure and Normalised Errors

MCSs are used to compare the performance of different QDGL parameter configurations. Table 2 shows the range in which the model parameters are randomised upon initialisation. Once the projectile and target are initialised with randomised speeds and positions, the simulation runs for a fixed time of 50s unless otherwise specified.

Table 2.

Range of values for initialisation of Monte Carlo simulation.

Parameter	Initialisation range
x₀, y₀	[ − 100, 100]
(u₀, v₀)	(0, 0)
ϕ ₀	[0, 2π]
xT₀, yT₀	[ − 10, 10]
(uT₀, vT₀)	(0, 0)

Multiple MCSs are run in a ‘batch’. Figure 6 shows the average system responses for four batches of 10⁴ MCSs. The difference in the response is caused by the stochastic nature of the MCS initialisation. Increasing the batch size results in an average system response which is sampled over a large number of simulations, providing more consistency between batches as well as being more representative of the true system behaviour. Obviously, the maximum discrepancy between multiple batches is inversely proportional to the batch size.

Figure 6.

Response discrepancy between multiple batches of 10⁴ MCSs.

Figure 7 shows the instantaneous maximum error for varying batch sizes over time. For a given batch size, four separate batches are run. The maximum discrepancy between any of the four batches at any specific instances is plotted, for duration of the simulations. This process is repeated for four different batch sizes.

Figure 7.

Instantaneous maximum variation of the normalised error of different batch sizes.

All batch sizes in the figure start with the same zero error, since the error curves which are being averaged all begin at the maximum error. In addition, the same QDGL parameters are used in every simulation in a given batch, so the identical terminal error reduction capabilities produce the same steady state error in all cases. This results in the maximum error difference being small at the beginning and end of every batch. Table 3 shows the peak error difference, total integral error and computation time for each batch size shown in Figure 7. A batch size of 10⁴ was chosen as the optimal trade-off between peak difference and computation time. While a batch size of 10³ produced an integral error of the same magnitude as 10⁴, the peak difference was unacceptably high.

Table 3.

Error characteristics and computation time for varying Monte Carlo simulation batch sizes.

MCSs per batch	Peak difference %	Total envelope error	Computation time
10²	15.56	2.95 × 10⁻¹	13 min
10³	4.70	7.16 × 10⁻²	2 h
10⁴	1.43	2.90 × 10⁻²	22 h
10⁵	0.42	6.88 × 10⁻³	9 Days

Each MCS produces a system response like that shown in Figure 8. The system response of each MCS is normalised against the initial distance error, which was in turn randomised at the beginning of each simulation. The instantaneous error, ϵ_t, is given by

ϵ_{t} = \frac{d_{t}}{d_{0}}

(18)

Figure 8.

System response showing different quasi-dynamic guidance law regimes.

The integral of the normalised system response is thus given by

ϵ = \int_{0}^{τ} \frac{ϵ_{t}}{ϵ_{0}} . d t

(19)

where є₀ is the initial error. By integrating the whole system response in this way, the system response for different QDGL configurations can be numerically compared. The size of the steady state error amplitude is small compared to the initial and transient error amplitudes. Thus, є_m is representative only of the initial convergence and transient errors. However, the error amplitude of overshoots are not negligible compared to transient error amplitudes; thus, overshoots will be detected by the integration. Likewise, the ‘batch average integral’ error,

\hat{ϵ}

, can be computed for a batch size M as

\hat{ϵ} = \sum_{m = 1}^{M} \frac{ϵ_{m}}{M}

(20)

where є_m is the normalised integral system response of MCS m. This value allows a direct and meaningful numerical comparison of batches with different QDGL parameter configurations.

Figure 8 shows an example system response with annotations showing the different regimes governed by equation (17). V_PT(d) = V_lim(d) − ξ is the closing regime, V_PT(d) = V_k is the linear regime and V_PT(d) = 0 is the stationary regime. In addition, the aspects of the steady state error are the transient, the steady state amplitude and the oscillation amplitude. The transient is also the motion of the projectile in the linear regime; a faster V_k provides a faster transient speed. The steady state amplitude can be thought of as the distance d of the ‘orbital centre’ from the target. The oscillation amplitude, and in addition the oscillation frequency, is governed by F_c and ω_0,1, in which neither can be affected by modifying V_k or ξ. These oscillations are caused by the holding orbit described in equation (12); if the locus of the assumed circular orbit coincides with the target, then there will be no oscillations. However, if any small perturbation offsets the orbit loci or the orbit isn’t perfectly circular, the amplitude of the steady state error and steady state oscillations will increase. If the steady state oscillation continues periodically, then the orbit is stable, and the projectile is remaining in the stationary regime. If the steady state has another lower frequency oscillation, the orbit is not stable. This is caused by the linear regime velocity V_k being too great for the current d₂, causing the projectile to pass straight through the stationary regime. N.B. The regimes are governed by d and so the vertical lines representing them on Figure 8 should be horizontal; this was intentional, to aid interpretation.

Figure 9 illustrates how V_lim is modified in the QDGL to account for the system response lag. V_lim acts as the reference signal which the QDGL then modifies using ξ. The closing speed function V_PT(d) is then passed to the autopilot, which calculates the change in speed per revolution demanded from the projectile’s actuator mechanism V_Dem. The V_Dem can then be passed through a chosen controller H, such as a PID. The V_Dem, modified under the action of controller H, is then passed to the actuator mechanism which saturates the signal such that V_Dem ∈ [0, ΔV_max]. N.B. V_max is an absolute limit of the system resulting from the maximum bias angle of one half roll revolution, and it is not a characteristic of any physical actuator hardware. If V_Dem > V_Max the system will saturate and will only be able to deliver V_Max. As such, any H > 1 when the autopilot is already demanding V_PT = V_max will have no effect on the system whatsoever.

Figure 9.

Modification of V_lim for speed controller with velocity feedback.

The sensor block shown in Figure 9 represents the sensor dynamics of the system which are assumed in the idyllic system to be perfect and instantaneous. For practical implementation, conventional on-board sensors or an external detection system could be used to measure the roll angle, similar to that described in Ref.18. If volumetric restrictions permit, on-board image sensing hardware could be used.³⁷ Linear ballistic theory could also be used to estimate the change in roll rate of the projectile along the trajectory beyond what is known from the projectile launch to supplement information from the sensors.

Parametric investigation

With the QDGL fully presented, the system parameters V_k, ξ, d₁, and d₂ can be investigated, with extreme variation in the parameters being shown to demonstrate their functionality and its impact on the system response. A more high-fidelity optimisation is conducted in the Results section. Throughout this section values such as ‘zero’, ‘low’, ‘medium’ and ‘high’ are used when discussing the system response. This is intended to illustrate the change in behavioural differences across the full range of suitable parametric values. However, the exact numeric values are given in the corresponding figure caption.

Figure 10 shows how modifying ξ to extreme values affects the system response. When ξ = 0, the system response exceeds V_lim which was the case in Figure 5 leading to a large overshoot and oscillatory motion where the projectile closes with a speed which is too high. Any negative value of ξ yields a steady state error that is unacceptably high; thus, it is not shown on the figure. For low values, 0 < ξ < 1, the error is reduced quickly with a low amplitude steady state error and no higher order oscillations. High values, ξ ≥ 1, produce a large amplitude steady state error.

Figure 10.

How extreme values of ξ affect the normalised error (zero: ξ = 0, low: ξ = 0.7, high: ξ = 2).

Figure 11 shows three example trajectories which correspond to the extreme values of ξ from Figure 10. The high steady state error for large ξ is represented by the projectile being brought to rest too far from the target. The high amplitude decaying oscillations in the normal error for ξ = 0 is visible as the projectile overshoots the target by a large margin. This is caused by too high a speed being demanded and the projectile is unable to reduce its speed in time due to the system lag. For low values of ξ, the projectile is brought sufficiently close to the target to initiate a regime change.

Figure 11.

How extreme values of ξ affect the picture plane trajectory (zero: ξ = 0, low: ξ = 0.7, high: ξ = 2).

Figure 12 shows the effect of modifying the linear regime constant velocity V_k. All speeds initially follow the same error reduction path from t = 0 to t = 12. This is the range governed by V_PT(d) = V_lim(d) − ξ and thus modifying V_k has no effect. If V_k is sufficiently small, as the projectile transitions from the linear velocity to the stationary regime, the orbital centre is brought to rest very close to the boundary of the stationary regime, d₂. This results in a higher amplitude steady state error than if the velocity was high enough to reduce the distance to d → 0 before it was brought to rest. This is apparent from the figure, as an increasing value of V_k results in a lower amplitude of the steady state error up to the point that V_k is too high resulting in an overshoot. The optimal V_k is a trade-off with d₂ to deliver the orbital centre sufficiently close to the target before switching to the stationary regime.

Figure 12.

How varying V_k affects the transient and steady state error (zero: V_k = 0, low: V_k = 0.01, medium V_k = 0.2, high: V_k = 1).

In this case, for V_k = 0, the linear regime vanishes, merging with the stationary regime, that is, V_PT(d|d ∈ [0, d₂]) = V_PT(d|d ∈ [d₂, d₁]) = 0. The result of this is that the projectile enters the stationary regime at a distance d₁ and this is apparent from the figure, with a steady but large magnitude steady state error. For low values of V_k, the transient is very slow, but the amplitude of the steady state error is small. For medium values of V_k, the most desirable system behaviour can be observed. There is a very quick transient period followed by a low steady state error amplitude. For high values of V_k, there is an unstable switch between the linear and stationary regimes, caused by a sufficiently high overshoot to exceed d₂.

Figure 13 shows the effect of modifying the boundaries of the linear regime, d₁ and d₁. While d₂ is already arbitrary and selected based on the chosen level of accuracy of the system, d₁ is calculated from d₁| (V_PT(d₁) − ξ = 0) (equation (15)), after a value of ξ has been selected. However, d₁ is varied manually here to illustrate the impact of linear regime size. If the boundaries are set to be equal, d₁ = d₂, then there is ‘no linear regime’, it is bypassed completely and the velocity is brought to relative rest immediately. This leads to a large steady error as was the case for low values of ξ and V_k. The steady state oscillation amplitude is the same as for any other case, since F_c, ω₀ and ω₁ are not being modified.

Figure 13.

How varying linear regime boundaries affects transient and the steady state error.

If d₂ is small then the projectile gets closer to the target before switching to the stationary regime, this is the ‘Late termination’. When d₂ is sufficiently small, it becomes significant compared to the distance that can be travelled by the projectile travelling at speed V_k during the time for one complete roll rotation. This results in an unstable steady state oscillation from overshooting, where the projectile continuously switches between the linear and stationary regime, which is indicated by the late regime termination on the figure. In an ‘early activation’, d₁ is higher than its true value would be when computed from equation (15). The projectile is brought to the linear regime speed V_k too early in time, at a point where V_lim(d) − ξ would have otherwise permitted a higher closing velocity, resulting in a transient period significantly longer than in the other cases. Desirable system behaviour is observed from the ‘ideal duration’ on the figure with a steady transient from the dynamic to the stationary regime followed by a stable steady state oscillation.

The steady state amplitude for early regime activation is lower than the ideal scenario, which is not expected since the lower bound of the linear regime is the same for both simulations. The stationary regime governs all d < d₂, which describes a circular area around the target of radius d₂. Since the projectile only calculates the bias points when ϕ = 0, the projectile will switch regimes at a different point depending on where the first bias calculation takes place within d₂. Small deviations in the trajectory can thus cause a discrepancy in steady state amplitude though the oscillation amplitude will remain the same in all cases. This discrepancy must be mitigated by averaging large number of simulations.

Results and discussion

This section discusses a genetic algorithm used to optimise the QDGL parameters, ξ and V_K, as well as the gains of a PID controller. As discussed in the previous section, the values of d₁ and d₂ cannot be optimised and are therefore not included in this section. The performance of the QDGL with the optimised ξ and V_K is then assessed for performance rejection capabilities and evaluated against stationary and moving targets.

System parameter optimisation using genetic algorithm

Genetic algorithms (GAs) are the proven method of optimising system parameters and PID controller gains.^38–41 A basic genetic algorithm is implemented using the parameters shown in Table 4. The GA was first tested on a drop-wave function of the form

f (x_{1}, x_{2}) = - \frac{1 + \cos (12 \sqrt{x_{1}^{2} + x_{2}^{2}})}{0.5 (x_{1}^{2} + x_{2}^{2}) + 2}

(21)

which is a non-convex, multi-modal, continuous function with

\min [f (x)] = - 1 when x = (0,0)

(22)

Table 4.

Range of values for initialisation of Monte Carlo simulation.

Parameter	Value
Generations G	200
Population size N	100
MC batch size M	10⁴

The GA converged to the same optimum value for 10 individual trial runs, to within 3sf, and verifying the GA can repeatedly converge to a known optimum solution in the given configuration.

Both the initialisation of the GA and the MCS procedure are stochastic in nature; hence, the optimum value found is not necessarily the ‘true’ optimum. However, the MCS procedure consistently represents the system response to within the desired degree-of-accuracy. In addition, the GA is of low dimensionality and is operating in a relatively low complexity space. This, in conjunction with the consistent performance of the GA when optimising the drop-wave function, reaffirms the optimum values provided by the GA are satisfactory for the scope of this paper.

Using this MCS batch procedure, each specimen of each generation within the GA can now be meaningfully compared such that an optimal solution may evolve. Algorithm 1 shows the order of operations for the GA. The fitness function, FIT, of the GA to be minimised is simply $\hat{ϵ}$ , since this represents the average error over time for many instances for a particular QDGL parameter candidate

FIT = \hat{ϵ}

(23)

The normalised integral error of the system response for the m^th MCS follows from equation (19) as

ϵ_{m} = \int_{0}^{τ} \frac{ϵ_{t}}{ϵ_{0}} . d t

(24)

Likewise, the mean normalised integral error of the system response for a Monte Carlo batch of size M for specimen n follows from equation (20) as

\hat{ϵ} = \sum_{m}^{M} \frac{ϵ_{m}}{M}

(25)

The Monte Carlo batch size for the proceeding is M = 10⁴, the justification for which was discussed in the previous section.

Algorithm 1

Execution of GA optimisation using the MCS procedure

1: Randomly initialise N specimens, ${(ξ, v_{K})}_{n}$

2: for Each generation g ∈ [1 : G] do

3: for Every Specimen in the population n ∈ [1 : N] do

4: Set QDGL parameters equal to specimen $(ξ, v_{K}) = {(ξ, v_{K})}_{n}$

5: Set $d_{1} = ξ^{2} / 2 \tilde{a}$ (From equation (16))

6: for For MCS m ∈ [1 : M] do

7: Run MCS m with random initial conditions

8: Compute ϵ_m

9: end for

10: Compute specimen Fitness: ${FIT}_{n} = \hat{ϵ_{n}}$

11: end for

12: Rank specimens in order of fitness and select candidates for reproduction

13: Create offspring from candidates and cull resulting population to size

14: Mutate population, then reduce mutation factor

15: end for

To reduce computation time, a preliminary search is conducted to inform the scope of the GA. Figure 14 shows a $\hat{ϵ}$ surface for the joint variation of ξ ∈ [0.2, one] and V_k ∈ [0, 1.5]. Low regions on the figure correspond to a low $\hat{ϵ}$ and the lowest point is with certainty, bounded by ξ ∈ [0.3, 0.8] and V_k ∈ [0, 0.5]. High values of ξ or V_K lead to a large $\hat{ϵ}$ . Large values of ξ produce a saturated response regardless of the value of V_k. In either case, V_PT(d) > V_lim(d) and there is a substantial overshoot, confirming what was seen in Figures 10 and 12. In addition, low values of V_k have a small rise in $\hat{ϵ}$ . This is not due to an overshoot, rather a slow transient results in a residual error which is not present for higher values of V_k which reduces the error quicker. In addition, this surface is less complex than the drop-wave function in which the GA was evaluated on during the preliminary tests, indicating that the GA will with high probability converge to the true solution.

Figures 15 and 16 show the convergence rate for the GA operating under the boundary conditions for ξ and V_k.

The optimal QDGL parameters from the GA were found to be ξ = 0.54 and V_k = 0.18.

Figure 14.

System response error $\hat{ϵ}$ as a function of quasi-dynamic guidance law parameters ξ and V_k.

Figure 15.

Convergence of 2D GA to optimise quasi-dynamic guidance law parameters.

Figure 16.

Generation distribution during initial convergence of GA for QDGL parameters.

PID controller gain optimisation

Figure 9 included a block H which represents the chosen controller for V_Dem. A PID controller is used to investigate the proportional, integral and differential aspects of V_Dem during feedback. The use may reveal system behaviour that was not otherwise obvious from the previously discussed framework and highlight any weaknesses of the QDGL approach. During the simulation, V_Dem is decomposed in the YZ plane to the YZ earth axis giving ${[u_{Dem} v_{Dem}]}^{T}$ . This vector is then modified by the chosen controller and exported to the actuator block which calculates the bias angles using the method discussed in the previous section. This means that one controller can act on both channels simultaneously or two controllers could act separately on each channel. Since the environment is simplified to a point of planar symmetry, one controller is chosen to act on both channels simultaneously. If external forces are introduced into the environment, then two controllers would be better suited to account for non-symmetrical forces which are biased to one direction, such as gravity. Using two controllers may require re-normalising the signal so any gain applied to one channel is not lost during saturation by the actuation mechanism.

Figure 17 shows a preliminary search of independently varying the PID controller gains, corollary to Figure 14. The gains are initially held at ${[k_{P} k_{I} k_{D}]}^{T} = {[\begin{matrix} 1 & 0 & 0 \end{matrix}]}^{T}$ , which emulates the absence of a PID controller. Each gain is individually swept over the given range while the other gains remain fixed, and a MCS batch is run to determine the corresponding $\hat{ϵ}$ for each specific gain configuration. As with the optimisation of the QDGL parameters, the intention is to reduce the scope of the GA optimisation and, therefore, computation time. The ranges bounding the optimum solution were found to be k_P ∈ [0.6, 1.4], k_I ∈ [0, 0.1] and k_D ∈ [0, 0.5]. The GA is then modified from algorithm 1, such that each specimen is now ${(k_{P} k_{I} k_{D})}_{n} \forall n \in N$ .

Figure 17.

Contribution of independent k_P, k_I and k_D to the system error.

Figures 18 and 19 show the convergence rate of the 3D-adapted GA. The optimal configuration of PID gains was found to be ${[\begin{matrix} k_{P} & k_{I} & k_{D} \end{matrix}]}^{T} = {[\begin{matrix} 1 & 0 & 0 \end{matrix}]}^{T}$ , indicating that tuning the QDGL parameters is a sufficient and complete optimisation for the described system. A PID controller will not provide a performance increase beyond the described tuning of the QDGL parameters, ξ and V_k, since there is no significant behavioural dependence on the state derivative or state-time integral. It is of note that the GA converged to a local minima much quicker than when optimising for ξ and V_k. This is likely due to the adverse effect any k_I, k_D > 0 has on the system response, which quickly coerces the evolution, effectively reducing the search to a 1D GA. The minimum $\hat{ϵ}$ achieved during the optimisation of the PID controller gains was equal to the minimum $\hat{ϵ}$ during the optimisation of the QDGL, to 3sf. This is within the expected MCS batch error.

Figure 18.

Convergence of 3D GA to optimise PID controller gains.

Figure 19.

Generation distribution during initial convergence of GA for PID controller gains.

Disturbance rejection and system performance

The optimised QDGL is now tested for disturbance rejection capabilities and performance against different target profiles. Figures 20 and 21 show how the projectile responds to different disturbances. In each case, the projectile and target are initialised at a specified distance, the target closes the distance under normal operation and is then allowed to remain in steady state for a sufficient time until the chosen disturbance is applied and synchronised at 50s.

Figure 20.

System response to various target distance disturbances.

Figure 21.

System response to various target velocity disturbances (low: V = 0.01 ms⁻¹, medium V = 0.1 ms⁻¹, high: V = 1 ms⁻¹).

Figure 20 shows disturbance displacements, where at 50s, the target coordinates are set to be at a magnitude of 0.5x, 1x and 2x that of the initial displacement, as indicated by the figure. For all magnitudes of displacement, the error change is discontinuous. The initial correction response is similar for all due to the demand of the velocity autopilot saturating the control mechanism. Once the projectile has slowed sufficiently, it enters the linear regime at the same point in each case, d ≈ 10⁻², since the regime switching is governed by a certain distance. There is a small discrepancy between linear regime switching for the disturbances and the initial reference signal. The reference signal enters slightly later at a lower distance. This is likely caused by the projectile crossing the regime threshold d₂ with more of the roll rotation left to complete, meaning it will travel longer before the speed is corrected again. The similarity in response is due to the fact that in all cases, the initial relative velocity between the projectile and target is zero, and thus, the system will respond as if the simulation has just been initialised at different distances.

Figure 21 shows velocity displacements, where at 50s, the relative velocity of the target is instantaneously changed to a low, medium and high respective speed, radially away from the projectile. The positions of the projectile and target are not changed; they are then free to dynamically evolve. When the velocity disturbance is low, within what the actuator mechanism is capable of correcting in one roll rotation, the disturbance is corrected quickly. With a medium disturbance, beyond the correction of one bias manoeuvre, the system takes longer to recover. Since this is a velocity disturbance, the maximum error increase is not instantaneous, rather it coincides with the instant where the target is no longer moving away from the projectile and the relative speed is zero. From this point, the closing of the projectile is similar to the distance disturbances. This is the same for the high velocity disturbance, except that the rate of reduction of error divergence takes longer to correct.

Figure 22 shows the average normalised error for a 10⁴ MCS batch against both stationary and moving targets. The response against stationary targets is the same as previously in this section. Against moving targets, however, the error initially increases a small amount before decreasing in a manner similar to the response against static targets. This initial increase is due to the random chance of the projectile being initialised with speeds directed away from the target, and then having to correct this dispersive motion before beginning the correction procedure.

Figure 22.

Performance of the quasi-dynamic guidance law against target profiles.

Conclusions

A 7 DOF dynamic model for a dual-spin projectile (DSP) is presented and implemented in computational simulations. A novel projectile design is presented along with the unconventional control method of asymmetric roll-rate biases. The quasi-dynamic guidance law (QDGL) is developed and a parametric study is conducted which shows how modifying QDGL parameters affect the system response. A Monte Carlo simulation (MCS) procedure is described, which is used to meaningfully compare the system response for different parameter configurations. A genetic algorithm which utilises the MCS procedure is then used to optimise the QDGL parameters and PID controller gains. The disturbance rejection capabilities of the optimised QDGL are tested, as well as the effectiveness against both static and dynamic targets. In all cases, the QDGL is able to reduce the distance error to a satisfactory level.

Full dynamic coupling and aerodynamic disturbances have not been considered in this paper but will be the topic of further investigations. In addition, the use of this QDGL should be explored with a dynamic or time dependant control force and arbitrarily complex roll rate switching profiles.

Footnotes

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This paper was sponsored by EPSRC ICASE Grant reference 1700 064 and BAE Systems

Data availability

The data that support the findings of this study are available from the corresponding author, James N., upon reasonable request.

ORCID iDs

James Norris

John Economou

References

Starstreak

. High velocity missile (hvm), 2016. Available at: https://www.thalesgroup.com/sites/default/files/database/document/2018-11/STARStreak_02_16SM.pdf

Global Security . XM1156 Precision Guidance Kit (PGK), 2011. Available at: https://www.globalsecurity.org/military/systems/munitions/m1156.htm

Silver

BAE

. Bullet – precision guidance kit. Available at: https://www.baesystems.com/en/product/silver-bullet–precision-guidance-kit 2015

Siouris

. Missile guidance and control systems. Heidelberg, Germany: Springer, 2004.

Wang

. Three–dimensional adaptive sliding mode guidance law for missile with autopilot lag and actuator fault. Int J Control Autom Syst 2019; 17(6): 1369–1377. DOI: 10.1007/s12555-017-0731-9

Geng

Zhang

Sun

. Adaptive back–stepping sliding mode guidance laws with autopilot dynamics and acceleration saturation consideration. Proc Inst Mech Eng G: J Aerospace Eng 2019; 233(13): 4853–4863. DOI: 10.1177/0954410019835728

Cui

, et al. Sliding mode control–based guidance law with damping term and angle constraint. Chin Control Conf 2019; 2019-July(61374160): 4083–4088. DOI: 10.23919/ChiCC.2019.8865893

Tan

. A novel sliding mode guidance law with impact angle and time constraints, 2019. 5896–5901.

Sharma

Ratnoo

. A bearings–only trajectory shaping guidance law with look-angle constraint. IEEE Trans Aerosp Electron Syst 2019; 55(6): 3303–3315. DOI: 10.1109/TAES.2019.2906090

10.

Lee

Kim

. Capturability of impact–angle control composite guidance law considering field–of–view limit. IEEE Trans Aerosp Electron Syst 2020; 56(2): 1077–1093. DOI: 10.1109/TAES.2019.2925485

11.

Wang

Dong

Wang

, et al. Nonlinear suboptimal guidance law with impact angle constraint: an sdre–based approach. IEEE Trans Aerosp Electron Syst 2020; 56(X): 4831–4840. DOI: 10.1109/TAES.2020.3003105

12.

Qiang

Jun

Yonghai

, et al. The research of terminal optimal guidance law of MANEUVERING vehicle with multiple constrains. In: IEEE 10th international conference on mechanical and aerospace engineering, Brussels, Belgium, 2274-2578 July 2019, pp. –. DOI: 10.1109/ICMAE.2019.8880990.

13.

Chwa

. Robust nonlinear disturbance observer based adaptive guidance law against uncertainties in missile dynamics and target maneuver. IEEE Trans Aerosp Electron Syst 2018; 54(4): 1739–1749. DOI: 10.1109/TAES.2018.2801392

14.

Alrawi

AAA

Graovac

Ahmad

, et al. Modified guidance law based on a sliding mode controller for a missile guidance system. Teh Vjesn 2016; 23(3): 639–646. DOI: 10.17559/TV-20130625173047

15.

Ratnoo

Shima

. Line–of–sight interceptor guidance for defending an aircraft. J Guidance, Control Dyn 2011; 34(2): 522–532. DOI: 10.2514/1.50572

16.

Yamasaki

Balakrishnan

Takano

. Modified command to line–of–sight intercept guidance for aircraft defense. J Guidance Control Dyn 2013; 36(3): 898–902. DOI: 10.2514/1.58566

17.

Zhe

Zhengjie

Ningjun

. Virtual target based optimal guidance law with impact angle and burst height constrains. Proc of ICMIC. In: 2015 7th international conference on modelling. identification and control, Sousse, Tunisia, 18-20 December 2015, pp. 1–6. DOI: 10.1109/ICMIC.2015.7409437.

18.

Khruslov

Feofilov

Goryachev

, et al. Missile control in the polar coordinate system using a scalar radius. Gyroscopy Navig 2015; 6(1): 66–72. DOI: 10.1134/S207510871501006X

19.

Penev

. An expanded two–dimensional proportional–derivative command to line–of–sight guidance law. Gyroscopy Navig 2018; 9(4): 344–351. DOI: 10.1134/S2075108718040132

20.

Zhang

Yao

Zheng

. Impact point prediction guidance based on iterative process for dual–spin projectile with fixed canards. Chin J Aeronautics 2019; 32(8): 1967–1981. DOI: 10.1016/j.cja.2019.06.002

21.

Guo

Song

Wang

, et al. Guidance law design for a class of dual–spin mortars. Int J Aerospace Eng 2015; 2015: 1–12. DOI: 10.1155/2015/952076

22.

Ratnoo

Ghose

. Collision–geometry–based pulsed guidance law for exoatmospheric interception. J Guidance Control Dyn 2009; 32(2): 669–675. DOI: 10.2514/1.37863

23.

Zarchan

. Tactical and strategic missile guidance. 6th ed. American Institute of Aeronautics and Astronautics Inc., 2012.

24.

Taur

Chern

. Passive ranging for dog–fight air–to–air IR missiles. In: 1999 guidance, navigation, and control conference and exhibit, 1999, pp. 1737–1751. DOI: 10.2514/6.1999-4289

25.

Song

Tae Yoon

. Practical guidance for homing missiles with bearings–only measurements. IEEE Trans Aerosp Electron Syst 1996; 32(1): 434–443. DOI: 10.1109/7.481284

26.

Theodoulis

Gassmann

Wernert

, et al. Guidance and control design for a class of spin-stabilized fin-controlled projectiles. J Guidance, Control Dyn 2013; 36(2): 517–531. DOI: 10.2514/1.56520

27.

Lee

Jun

. Autopilot design for dual–spin projectile based on PI and feedback linearization control. IEEE conference on control applications,

IEEE, Juan Les

Antibes, France, 82084-102088 October 2014, pp. –. DOI: 10.1109/CCA.2014.6981610

28.

Seve

Theodoulis

Wernert

, et al. Flight dynamics modeling of dual–spin guided projectiles. IEEE Trans Aerosp Electron Syst 2017; 53(4): 1625–1641. DOI: 10.1109/TAES.2017.2667820

29.

Fan

. Research on instability boundaries of control force for trajectory correction projectiles. Math Probl Eng 2019; 2019: 1–12. DOI: 10.1155/2019/6362835

30.

McCoy

. Modern exterior ballistics: the launch and flight dynamics of symmetric projectiles. Pennsylvania, USA: Schiffer Publishing, 1999.

31.

Lloyd

Brown

. Instability of spinning projectiles during terminal guidance. J Guidance Control 1979; 2(1): 65–70. DOI: 10.2514/3.55833

32.

Robinson

JWC

Berefelt

. On guidance and control for guided artillery projectiles, part 1: general considerations, 2011. (October). Available at: http://foi.se/ReportFiles/foir_3291.pdf.

33.

Norris

Hameed

Economou

, et al. A review of dual–spin projectile stability. Defence Technol 2020; 16: 1–9. DOI: 10.1016/j.dt.2019.06.003

34.

Fairfax

Vasile

Strohm

, et al. Trajectory shaping for quasi-equilibrium glide in guided munitions 2020, 2020, pp. 1–12. DOI: 10.2514/6.2020-0021

35.

Wang

Cheng

, et al. Influence of yawing force frequency on angular motion and ballistic characteristics of a dual–spin projectile. Defence Technol 2016; 12(2): 124–128. DOI: 10.1016/j.dt.2015.12.007

36.

Wang

. Influence of control strategy on stability of dual–pin projectiles with fixed canards. Defence Technol 2018; 14(6): 709–719. DOI: 10.1016/j.dt.2018.04.014

37.

Fresconi

Rogers

. Flight control of a small–diameter spin–stabilized projectile using imager feedback. J Guidance Control Dyn 2015; 38(2): 181–191. DOI: 10.2514/1.G000815

38.

Lin

Bao

Braasch

, et al.

Genetic algorithm based parameter optimization of a fuzzy logic controller. 1999 guidance, navigation, and control conference

and exhibit 1999, Portland, OR, U.S.A, 09 August 1999 - 11 August 1999, 1117–1122. DOI: 10.2514/6.1999-4152.

39.

Eun

Bang

. Cooperative task assignment and path planning of multiple UAVs using genetic algorithm. In: Collection Tech Pap – 2007 AIAA InfoTech Aerospace Conference, Rohnert Park, California, 3, May 2007: 2680–2697. DOI: 10.2514/6.2007-2982.

40.

Hong

Pey

Soh

. Structural parameter identification using robust genetic algorithm optomization method, 2004. DOI: 10.2514/6.2004-2015.

41.

Kramer

. Studies in Computational Intelligence 679 Genetic Algorithm Essentials. ISBN 9783319521558. 2017.