AIAA 2013-5114

Perturbation Based Guidance for a Generic 2D Course Correcting Fuze John W.C. Robinson* and Peter Str¨omb¨ack†

Downloaded by John Robinson on September 5, 2013 | http://arc.aiaa.org | DOI: 10.2514/6.2013-5114

Dept. Information and Aerosystems, Swedish Defence Research Institute (FOI), Stockholm, Sweden

A guidance strategy for generic 2D course correcting fuzes with canard control is developed using a perturbation approach. The starting point is a nominal trajectory which is computed before flight. During flight, only simple correction terms need to be calculated in order to form the guidance commands. The resulting strategy is of essentially predictive character but can be used with shorter prediction horizon than time-to-go, thereby effectively implementing a path following strategy. Simulations are provided to illustrate the effectiveness of the approach for a detailed model of a standard 155mm artillery shell in varying firing scenarios. It is shown that dispersion figures in the order of meters are possible to obtain with realistic values for errors in muzzle velocity, gun aiming and available control authority in the canard control system.

I.

Introduction

The extended range that indirect fire systems offer over direct fire systems comes at a price in the form of increased dispersion. The range achievable with a typical 155mm artillery system is well over 20km, even without using base bleed for drag reduction or rocket assist to boost range. At these ranges the dispersion in impact point of standard unguided rounds can be several times the kill radius (about 50m) of a typical high explosive fragmentation type shell, at least in the downrange direction. With base bleed or rocket assist the dispersion becomes even greater. Therefore, it has been natural to try to bring down the dispersion by deployment of specially designed guided projectiles or adding guidance and control functions to traditional unguided rounds.1 The former line of development has lead to precision guided projectiles such as the Raytheon/BAE Systems M982 Excalibur (GPS navigation) and the KBP KM-2 Krasnopol (laser target designation). These types of projectiles, although offering dispersions in the order of meters (in their latest versions), are very costly compared to standard unguided rounds. Moreover, the precision offered by these advanced projectiles is frequently higher than required by the mission (e.g. a stationary target, extending over tens of meters). A compromise solution that can address both these problems is provided by course correcting fuzes. Course correcting fuzes can be mounted onto standard, existing rounds and have the capability of reducing the impact dispersion to acceptable values at only a moderate increase in cost. I.A.

Controlling dispersion with course correcting fuzes

The major contributing factors to dispersion of unguided artillery projectiles are deviations from nominal values in muzzle velocity and meteorological conditions (wind, air density, temperature), aiming errors of the gun and deviation in projectile properties. Downrange dispersion is largest and can be many times larger than the crossrange dispersion, for long ranges. Therefore, some of the existing course correcting fuzes employ mainly downrange correction, using deployable drag brakes as the main means of control.2–4 When controlling the downrange errors is not enough, drag brakes are often combined with some means of crossrange correction, such as spin brakes which reduce the yaw of repose, hence the crossrange drift. Still, spin brakes yield only small crossrange correctional capabilities.5 For this reason, 2D (downrangecrossrange) correction techniques based on canard control, which is a natural next step in complexity, have been introduced in course correcting fuzes. * Deputy † Senior

Research Director, AIAA member, [email protected] Scientist, AIAA member, [email protected]

1 of 26 American Institute of Aeronautics and Astronautics Copyright © 2013 by John W.C. Robinson, Peter Strömbäck. Published by the American Institute of Aeronautics and Astronautics, Inc., with permission.

The 2D canard control concept is not new, however. A design based on a combination of movable and fixed canards was proposeda and analyzed in the 1970’s by Regan and Smith.7 In this design, one pair of canards are set at a (small) fixed cant angle to make the fuze rotate slowly (controlled by a spin brake) in the opposite direction of the shell while the remaining pair of canards can be deflected to produce a net normal force in any desired direction. Variants of this idea have subsequently been implemented in existing course correcting fuzes.8 We shall here consider a generic variant of the 2D canard control concept, where it is supposed that the normal force exerted by the canards at the nose can be directed at any direction at any time and continuously adjusted,9 thus neglecting to model how this is actually accomplished.

Downloaded by John Robinson on September 5, 2013 | http://arc.aiaa.org | DOI: 10.2514/6.2013-5114

I.B.

Guidance in course correcting fuzes

Even though 2D canard control can offer a high degree of deflectability to a nominal ballistic trajectory, there are natural drawbacks and limitations with canard control of spinning projectiles. One drawback is that a step in normal force applied some distance forward of the center of mass (CoM) on a spinning projectile yields an “out-of-phase” swerve response whereas a force applied near, or aft of, the CoM does not.10 In fact, for spinning projectiles it is preferable to apply the normal force at the base of the projectile since it yields the best control authority.10, 11 Furthermore, a too large normal force near the tip of the projectile can easily lead to lossb of stability.13, 14 A consequence of this is that essentially only (small) course correction (as opposed to larger trajectory shaping) can be obtained with spinning airframes and canard control, a fact which has lead to renewed interest15 in fin stabilized airframes also for projectiles launched from gun barrels with twist.c This also raises the question of how large system or model errors can be corrected with course correcting fuzes based on canard control, a question which we shall address in this work. The three main types of guidance most often considered in guided munition contexts are trajectory shaping, trajectory prediction and path following.17 When only small corrections to a trajectory are possible, only the latter two are applicable. In predictive guidance a model of the projectile and the environment is used in each guidance update instant to compute a sequence of control actions which would yield a flight path leading to the target, given the current state.18 If no errors exist in the model, the state is perfectly known and actuators are ideal, a single guidance computation and following control action is sufficient to give the projectile a flight path leading to the target. However, in reality at least a few such computations with accompanying corrections of the path are needed during the flight, hence the inherently open loop predictive guidance strategy will in practice lead to a closed loop one.d Predictive guidance can also be applied when the prediction horizon is shorter than the estimated time-to-go to the target and the desired terminal state is instead a point on the flight path. This can be thought of as a form of path following guidance. In path following guidance, often called trajectory following guidance,e the current position of the projectile is compared with the nominal flight path of the projectile and a control is applied in order to bring the projectile back to the nominal path in some near future (i.e. before it reaches the target). A drawback with path following guidance, however, is the fact that it only uses the position to correct the flight path. This means that each time a correction leads the projectile back to the nominal path it does not in general have right velocity to continue to follow the path but flies past it. Hence, further corrections are inevitable and higher control effort is required. I.B.1.

Related work

Both predictive guidance and path following guidance for spinning projectiles have received some interest in the literature. Gagnon and Lauzon17, 20 investigated path following guidance techniques applied to a standard 155mm artillery projectile. They studied three different types of actuation; drag brake, spin brake and roll decoupled fuze with four canards. It was found17 that drag brakes could effectively correct muzzle a An

early aerodynamic investigation of a 105mm projectile with movable canards can be found in Ref. 6. reason for loss of stability, in a real implementation with roll stabilized fuze, is that the axial moment induced between the fuze and the body can destabilize the projectile.12 c Hybrid concepts of spin and fin stabilized projectiles have also been considered, see Ref. 16 where the low drag of a spin stabilized body is combined with the end game maneuverability of a fin stabilized body with deployable canards for control. d For maneuvering targets the most common guidance method is proportional navigation19 which is partially predictive and uses current values of velocity and possibly acceleration to compute a control that successively decreases the miss predicted distance. e Here we distinguish between the path of a trajectory (the collection of all position points; the flight path) and the trajectory itself (the map from time to position-velocity pairs). b Another

2 of 26 American Institute of Aeronautics and Astronautics Copyright © 2013 by John W.C. Robinson, Peter Strömbäck. Published by the American Institute of Aeronautics and Astronautics, Inc., with permission.

Downloaded by John Robinson on September 5, 2013 | http://arc.aiaa.org | DOI: 10.2514/6.2013-5114

velocity errors but canards were superior for aiming errors and cross wind errors. Predictive guidance for direct fire spinning projectiles was investigated by Calise and El-Shirbiny.21 It was observed that the feedback path formed by the guidance loop can cause instability for this type of guidance due to excitation of the precession mode. Their prediction algorithm was based on a modified point mass model (MPM) in the vertical plane and on a second order polynomial approximation in the horizontal plane. Ollerenshaw and Costello18 studied predictive guidance for a direct fire, fin stabilized projectile with canards, in the setting of discrete time linear quadratic control. The impact point was predicted in-flight using so called projectile linear theory. It was found that the control could dramatically decrease the dispersion due to launch errors and that the discretization step of the flight path had little effect on the dispersion, as long as it was sufficiently small. The prediction horizon was also found to have a large effect on the results, with a longer horizon leading to reduced dispersion (at the cost of increased computational complexity). The stability and disturbance rejection for 2D canard controlled spinning projectiles can be improved by incorporation of a dynamics synthesizing autopilot and this has been investigated by Theodoulis and Wernert,22 and very recently by Theodoulis et al.23 They have shown that good tracking performance can be achieved by a controller of the proportional-integral type designed using linear parameter varying theory incorporating ℋ2 or ℋ∞ conditions. A dynamics synthesizing autopilot will inevitably have a bandwidth beyond that of the precession frequency, however, and thus requires fast actuators. The design of a controller for precession control using canard actuators that oscillate at the precession frequency is described in Ref. 24. I.B.2.

Contribution

Predictive guidance has the potential of being very energy efficient and requiring low control effort but it requires a good prediction model. A basic requirement for a course correcting fuze on a 155mm projectile can be to bring down the dispersion to about half the kill radius of a typical (25m) shell since this would mean that the target would be destroyed with high probability in a single shot. However, for long prediction horizons (≥ 25s, corresponding to more than half of the time of flight), only a full 6DOF model or a model based on modified projectile linear theory is in general capable25 of producing impact point prediction errors of this size. In the present work we propose a predictive guidance strategy that aims to overcome this obstacle. It is based on perturbation theory using a nominal trajectory that can be computed off-line, with arbitrarily high precision. Then, during flight, only low-complexity perturbation terms relating to this trajectory need to be computed in order to form the guidance command. The result is a guidance strategy that retains the low complexity of the path following strategies while at the same time offering the performance and energy economy of a predictive strategy. We show in particular that very simple corrections, corresponding to the ones needed for a perturbed vacuum trajectory, will provide very efficient dispersion reduction. The strategy is illustrated with simulations of a M107 projectile with a generic 2D course correcting fuze with canard control. I.C.

Notation

Notation is standard and we shall largely omit units. Points in R𝑛 × R𝑛 will sometimes be denoted as (𝑥, 𝑦) and sometimes as [𝑥𝑇 , 𝑦 𝑇 ]𝑇 where superscript 𝑇 denotes transpose. The norm ‖ · ‖ is always the 2-norm. I.D.

Outline

In the next section we introduce some notation and describe our trajectory model and its associated approximations. After this, the guidance strategy is described in section III. The simulations are given in section IV and some conclusions are offered in section V. Some calculations are deferred to appendix, such as the calculations of higher order approximations of perturbation terms and expressions for the autopilot used in the simulations.

II.

Trajectory Model

In order to accurately predict projectile trajectories, nonlinear models with six degrees of freedom (6DOF) and twelve state variables are normally required (for the rigid body motion of the shell). A good approximation to the 6DOF solution can however be obtained by the modified point mass trajectory model which has four degrees of freedom and seven state variables.26 The MPM model is based on the assumption that

3 of 26 American Institute of Aeronautics and Astronautics Copyright © 2013 by John W.C. Robinson, Peter Strömbäck. Published by the American Institute of Aeronautics and Astronautics, Inc., with permission.

the (uncontrolled) projectile at all times flies at equilibrium for the (vectorial) yaw of repose. This neglects the details of the oscillation of the projectile around the yaw of repose, i.e. the precession and nutation, and has differential equations only for the motion of the CoM and the spinning motion around the main axis. However, there is no condition on the spin rate imposed by the firing scenario, i.e. gun and target location, and the spin rate evolution is only weakly affected by small perturbations of the trajectory. Hence, the spin rate may effectively be considered as a function of time (for small perturbations around a nominal trajectory). With this simplification the MPM model can be thought of as being time varying but having only three degrees of freedom and six state variables. In the present work it is assumed that such models are sufficient to model the projectile trajectories for course corrected projectiles (this is further motivated below). The guidance equations developed are based on a generic form of the equations of motion with six states (position and velocity) which contains the time varying simplified MPM model with six states as a special case.

Downloaded by John Robinson on September 5, 2013 | http://arc.aiaa.org | DOI: 10.2514/6.2013-5114

II.A.

Firing scenario

The most natural frame of reference for guidance of projectiles is the fire control frame 𝐹 . This frame is a Cartesian right handed frame with the 𝑥-direction pointing downrange, the 𝑦-direction pointing upwards and the 𝑧-direction pointing crossrange to the right,f as seen when looking forward, in the direction of the barrel of the gun. The fire control frame is the (implicit) frame of reference for the discussion on trajectory modeling and the guidance equations developed below. We assume that the muzzle of the gun is located at a point with coordinates 𝑥𝐺 ∈ R3 referred to 𝐹 and the target (desired terminal point) is located at a (fixed) point in space given by 𝑥𝑇 ∈ R3 referred to 𝐹 (neither 𝑥𝐺 nor 𝑥𝑇 need represent a point on the Earth’s surface). We assume further that the target location relative to the gun is such that a firing solution exists, which we define next. II.B.

Firing solution

It will be important to emphasize dependence on initial values and therefore we shall adopt the formalism (mainly notation) of mathematical (vector field) flows to describe the projectile trajectories. II.B.1.

Firing flow

Given the position of the gun 𝑥𝐺 ∈ R3 , a firing flow Φ𝑥𝐺 is a (𝐶 1 , possibly time dependent) flow Φ𝑥𝐺 : ℳ × ℐ → ℳ defined on a manifold ℳ ⊆ R6 and for an open time interval ℐ containing 0, such that (𝑥𝐺 , 𝑣 0 ) ∈ ℳ for some initial muzzle velocity 𝑣 0 ∈ R3 . A firing flow thus represents the family of solutions to the differential equations used to describe the motion of the CoM when the model has six states and the initial conditions are varied. Given also the position of the target 𝑥𝑇 ∈ R3 , a firing flow Φ𝑥𝐺 is a firing solution flow ifg Φ𝑥𝐺 ((𝑥𝐺 , 𝑣 𝐺 ), 𝑡𝑓 ) = (𝑥𝑇 , 𝑣 𝑓 ) (1) for some 𝑣 𝐺 ∈ R3 and 𝑣 𝑓 ∈ R3 , 𝑡𝑓 ≥ 0, where 𝑡𝑓 is an unspecified (free) terminal time and 𝑣 𝑓 is the unspecified (free) terminal velocity. A firing solution is a triple (𝑥𝐺 , 𝑥𝑇 , 𝑣 𝐺 ) such that (1) holds for some terminal time 𝑡𝑓 and terminal velocity 𝑣 𝑓 . Nominal trajectory.

A (projectile, CoM) trajectory is a map 𝑡 ↦→ Φ𝑥𝐺 ((𝑥0 , 𝑣 0 ), 𝑡),

where Φ𝑥𝐺 is a firing flow, and a nominal trajectory is a firing solution flow Φ𝑥𝐺 evaluated for an initial velocity 𝑣 0 = 𝑣 𝐺 such that (1) holds. In other words, a firing solution flow Φ𝑥𝐺 represents a family of solution trajectories 𝑡 ↦→ Φ𝑥𝐺 ((𝑥, 𝑣), 𝑡), 𝑡 ∈ ℐ, (𝑥, 𝑣) ∈ ℳ, to the six-state equations of motion for the projectile, such that (at least) the particular solution starting at (𝑥𝐺 , 𝑣 𝐺 ) ∈ ℳ passes through the target location 𝑥𝑇 at time 𝑡 = 𝑡𝑓 . Since the flow Φ𝑥𝐺 is continuously f The exact orientation of the coordinate axes in 𝐹 relative to the gun can vary slightly. One choice is to have the 𝑥-axis exactly aligned with the gun, when the gun is pointing in exactly the right direction (i.e. no azimuth error). g Note that 𝑣 𝐺 in (1) is in general not unique; compare the the high and low elevation solutions for non maximal range fire. (Neither 𝑡𝑓 , 𝑣 𝑓 need be unique, in a general case, but we shall tacitly assume that they are.)

4 of 26 American Institute of Aeronautics and Astronautics Copyright © 2013 by John W.C. Robinson, Peter Strömbäck. Published by the American Institute of Aeronautics and Astronautics, Inc., with permission.

differentiable, the position-velocity pair (𝑥1 , 𝑣 1 ) for the projectile at some later time 𝑡1 > 0 is a continuously differentiable function of the initial condition values (𝑥0 , 𝑣 0 ). The firing flow formalism developed here is thus an abstract but simple way to describe generic trajectory models with three degrees of freedom and six states which emphasizes their “forward map” property, i.e. the dependence on initial values.h Moreover, the formalism extends naturally to the situation where there is an external control input 𝑢 present which acts as an extra parameter for the flow. In the following sections we shall use an explicit representation of such a controlled flow to develop the guidance equations. II.C.

State space representation and controlled dynamics

Downloaded by John Robinson on September 5, 2013 | http://arc.aiaa.org | DOI: 10.2514/6.2013-5114

The controlled flow we shall use to model the trajectories of a projectile (with course correcting fuze) is given by the solution to the ordinary differential equation {︃ 𝑥˙ = 𝑣 , 𝑡 ∈ ℐ, (2) 𝑣˙ = 𝑔 + 𝑎 + 𝑢 where 𝑔 ∈ R3 represents gravity (constant) and 𝑎 : ℳ × ℐ → R3 , which is 𝐶 1 (at least one time continuously differentiable in its arguments), represents the aerodynamic acceleration at free flight trim, i.e. yaw of repose (with control surfaces in their nominal setting, corresponding to 𝑢 = 0), and possibly Coriolis acceleration. The control actionsi are represented in open loop by the acceleration term 𝑢 : ℐ → R3 , which is 𝐶 1 . (In closed loop 𝑢 becomes a function of the state (𝑥, 𝑣) and possibly time.) In other words, we assume a control affine model. This means that the control acceleration 𝑢 in the model (2) does not only physically represent the control forces exerted by the nose mounted canards. Instead, 𝑢 represents the net resulting normal force resulting from the altered trim state (away from the yaw of repose) due to deflection of the canards.j Since 𝑎 is nonlinearly dependent on the position-velocity pair (𝑥, 𝑣) the system of differential equations in (2) is nonlinear but it has a unique solution (at least locally in time, for some open interval ℐ containing 0) for any initial condition (𝑥0 , 𝑣 0 ) ∈ ℳ since 𝑎 is 𝐶 1 . The family of solutions obtained by varying the initial conditions defines a flow Φ : ℳ × ℐ → ℳ. Remark: The model (2) is more general than it may first seem. Since (2) can contain any type of empirical corrections to a given nominal trajectory, it can in practice be made to arbitrarily well reproduce any trajectory (i.e. position-velocity history) coming from a more physically correct model, such as a 6DOF model. Hence, for the development of guidance equations, we may consider the real physical trajectory for the CoM (or a 6DOF model of it) to be a solution to an equation of the form (2). II.C.1.

Controlled firing flow

A firing (solution) flow with a control acceleration added, as for the solution to (2), will be called a controlled firing (solution) flow. In order to indicate the parametric dependence in Φ on the control function 𝑢( · ) we shall denote the solution to (2) starting at time 0 from the point (𝑥0 , 𝑣 0 ) and evaluated at time 𝑡 by Φ((𝑥0 , 𝑣 0 ), 𝑡; 𝑢). A solution obtained for 𝑢 = 0 (i.e. with actuators fixed in the zero position; an uncontrolled solution) will be denoted Φ((𝑥0 , 𝑣 0 ), 𝑡), i.e. we set Φ((𝑥0 , 𝑣 0 ), 𝑡) = Φ((𝑥0 , 𝑣 0 ), 𝑡; 0). In the applications below the initial point (𝑥0 , 𝑣 0 ) can represent any point along a firing flow trajectory (hence it need not signify the gun position and muzzle velocity). Therefore, the time 𝑡 = 0 can refer to any point along the trajectory which is considered as starting point for a certain part of the discussion. II.C.2.

Representation as linear perturbation of the vacuum trajectory

The solution trajectory 𝑡 ↦→ Φ((𝑥0 , 𝑣 0 ), 𝑡; 𝑢) to (2) is the same as what would be obtained if we substituted ¯ for 𝑎 + 𝑢 on the right in (2), where 𝑎 (︀ )︀ ¯ ((𝑥0 , 𝑣 0 ), 𝑡; 𝑢) = 𝑎 Φ((𝑥0 , 𝑣 0 ), 𝑡; 𝑢), 𝑡 + 𝑢. 𝑎 (3) h Being flows, the firing flows of course have a transitive property (see below) and also a “backwards map” property, i.e. dependence on terminal values. i Even though we shall eventually consider fuzes that only have 2D control actuation we let 𝑢 take values in R3 in the discussion here focusing on guidance. j This must be kept in mind when considering the equations for the autopilot in Appendix.

5 of 26 American Institute of Aeronautics and Astronautics Copyright © 2013 by John W.C. Robinson, Peter Strömbäck. Published by the American Institute of Aeronautics and Astronautics, Inc., with permission.

¯ which (That is, we first solve (2) as it stands and then use the solution to define a time varying quantity 𝑎 can be used instead of 𝑎 + 𝑢 when solving (2); the result (for fixed initial conditions) of the two approaches ¯ is (unlike 𝑎) only directly dependent on time, on the initial condition will be the same.) The function 𝑎 (𝑥0 , 𝑣 0 ) and on the choice of control function 𝑢( · ) (both in the open and closed loop cases). Furthermore, ¯ is (by 𝐶 1 dependence of solutions to ODEs with respect to initial conditions) 𝐶 1 in the initial the function 𝑎 condition (𝑥0 , 𝑣 0 ). Thus, we can represent the solution Φ to the system (2) by the (time invariant, vectorial) variation-of-parameters formula [︃ ]︃ ∫︁ 𝑡 𝑥0 ¯ ((𝑥0 , 𝑣 0 ), 𝑠; 𝑢)) 𝑑𝑠, 𝑡 ∈ ℐ, Φ((𝑥0 , 𝑣 0 ), 𝑡; 𝑢) = exp(𝑡𝐴) + exp((𝑡 − 𝑠)𝐴)𝐵(𝑔 + 𝑎 (4) 𝑣0 0 where

[︃

Downloaded by John Robinson on September 5, 2013 | http://arc.aiaa.org | DOI: 10.2514/6.2013-5114

𝐴=

03×3 03×3

𝐼 3×3 03×3

]︃

[︃ ,

𝐵=

03×3 𝐼 3×3

]︃ (5)

(i.e. forced double integrator dynamics). The formula (4) is a special case of a relation for perturbed nonlinear systems due to Alekseev27 (see also Ref. 28). We note here that the matrix 𝐴 is nilpotent, more precisely 𝐴2 = 06×6 , and therefore the series defining the matrix exponential contains only two terms, viz. [︃ ]︃ 𝐼 3×3 𝜉𝐼 3×3 exp(𝜉𝐴) = 𝐼 6×6 + 𝜉𝐴 = . 03×3 𝐼 3×3 It follows that ∫︁

∫︁ 𝑡 [︃

𝑡

exp((𝑡 − 𝑠)𝐴)𝐵𝑔 𝑑𝑠 = 0

0

(𝑡 − 𝑠)𝐼 3×3 𝐼 3×3

]︃

[︃ 𝑑𝑠 𝑔 =

𝑡2 2 𝐼 3×3

𝑡𝐼 3×3

(6)

]︃ 𝑔.

(7)

¯ obtained for 𝑢 = 0 will be denoted 𝑎 ¯ ((𝑥0 , 𝑣 0 ), 𝑡), i.e. we set The acceleration 𝑎 ¯ ((𝑥0 , 𝑣 0 ), 𝑡) = 𝑎 ¯ ((𝑥0 , 𝑣 0 ), 𝑡; 0). 𝑎

(8)

For zero air density (and nonrotating Earth) the term 𝑎 in (2) is identically zero and we obtain the basic ballistic vacuum trajectory as the solution (with 𝑢 = 0). II.D.

Perturbations of the nominal trajectory

The fact that an uncontrolled firing flow Φ, for instance represented by the solution (4) to (2) for 𝑢 = 0, is assumed to be continuously differentiable in its arguments means in particular that for ((𝑥0 + ∆𝑥, 𝑣 0 + ∆𝑣), 𝑡 + ∆𝑡) in a neighborhood of ((𝑥0 , 𝑣 0 ), 𝑡) ∈ ℳ × ℐ we have Φ((𝑥0 + ∆𝑥, 𝑣 0 + ∆𝑣), 𝑡 + ∆𝑡) = Φ((𝑥0 , 𝑣 0 ), 𝑡) (︀ )︀ 𝜕Φ((𝑥, 𝑣 0 ), 𝑡) ⃒⃒ 𝜕Φ((𝑥0 , 𝑣), 𝑡) ⃒⃒ 𝜕Φ((𝑥0 , 𝑣 0 ), 𝑠) ⃒⃒ + ∆𝑥 + ∆𝑣 + ∆𝑡 + 𝑜 ‖((∆𝑥, ∆𝑣), ∆𝑡)‖ , (9) 𝑥=𝑥 𝑣=𝑣 𝑠=𝑡 0 0 𝜕𝑥 𝜕𝑣 𝜕𝑠 where the multidimensional partial derivatives denote Jacobians and 𝑜 is (Landau’s) small o (ordo). Thus, for a sufficiently small perturbation ((∆𝑥, ∆𝑣), ∆𝑡) in initial position-velocity and terminal time, the change ∆Φ in the value of the flow (i.e. position-velocity) at the perturbed terminal time 𝑡 + ∆𝑡 due to the perturbation (∆𝑥, ∆𝑣) in initial values is well approximated by its first order (linear) approximation, i.e. by the three perturbation terms on the right in (9). The guidance principle we propose below is based on this simple property applied to (4). II.D.1.

Coriolis acceleration.

It is very easy to write down an exact expression for the Coriolis contribution to 𝑎. However, the effect is very small, and moreover linear in 𝑣, and therefore will essentially vanish in the perturbation expressions for the fire solution flow (since the net effect is of size 𝜔𝐸 ‖∆𝑣‖, where 𝜔𝐸 is the magnitude of Earth’s angular velocity about the north-south pole axis). Hence, we shall neglect the Coriolis contribution and therefore 𝑎 will henceforth represent only aerodynamic acceleration. 6 of 26 American Institute of Aeronautics and Astronautics Copyright © 2013 by John W.C. Robinson, Peter Strömbäck. Published by the American Institute of Aeronautics and Astronautics, Inc., with permission.

II.D.2.

Form of the perturbation terms

To study in more detail the effect of perturbations of the initial condition (𝑥0 , 𝑣 0 ) on the solution Φ in (4) to (2) in the uncontrolled case (𝑢 = 0) we can note that [︃ ]︃ ⃒ ∫︀ 𝑡 0 ),𝑠) ⃒ 𝐼 3×3 + 0 (𝑡 − 𝑠) 𝜕 𝑎¯ ((𝑥,𝑣 𝑑𝑠 𝜕Φ((𝑥, 𝑣 0 ), 𝑡) ⃒⃒ 𝜕𝑥 𝑥=𝑥 0 ∫︀ 𝑡 𝜕 𝑎¯ ((𝑥,𝑣0 ),𝑠) ⃒ = (10) 𝑥=𝑥0 ⃒ 𝜕𝑥 𝑑𝑠 𝜕𝑥 0 𝑥=𝑥0

and

Downloaded by John Robinson on September 5, 2013 | http://arc.aiaa.org | DOI: 10.2514/6.2013-5114

𝜕Φ((𝑥0 , 𝑣), 𝑡) ⃒⃒ = 𝑣=𝑣 0 𝜕𝑣

[︃

∫︀ 𝑡

]︃ ⃒ 0 ,𝑣),𝑠) ⃒ (𝑡 − 𝑠) 𝜕 𝑎¯ ((𝑥𝜕𝑣 𝑑𝑠 𝑣=𝑣 0 ⃒ ∫︀ 𝑡 0 ,𝑣),𝑠) ⃒ 𝐼 3×3 + 0 𝜕 𝑎¯ ((𝑥𝜕𝑣 𝑑𝑠 𝑣=𝑣 0

𝑡𝐼 3×3 +

0

(11)

(where differentiation under the integral sign is easy to justify). From this it follows in particular that the ¯ is quadratic in 𝑡, for small 𝑡. As a rough approximation (but effect of the aerodynamic acceleration term 𝑎 in practice surprisingly useful, see below) one can use the vacuum trajectory version of the expressions in ¯ in the integrals in (10) and (11) to zero. For more accurate (10) and (11), which is obtained by setting 𝑎 approximations, the contribution from the integrals must be taken into account. Since the integrands in the integrals in (10) and (11) are defined in terms of the nominal trajectory, the integrals can be computed offline and stored in tabulated form in the guidance computer. However, there might be situations where such precomputation and storage is not possible and therefore it is required to have also easily computed approximations for the integrals available. Moreover, there are alternative methods to direct computation (numerical differentiation combined with integration) for computing the approximations, as we show next. For the effect of a perturbation in the terminal time we get similarly in the uncontrolled case that [︃ ]︃ 𝜕Φ((𝑥0 , 𝑣 0 ), 𝑠) ⃒⃒ [Φ((𝑥0 , 𝑣 0 ), 𝑡)]4:6 = , (12) 𝑠=𝑡 𝜕𝑠 ¯ ((𝑥0 , 𝑣 0 ), 𝑡) 𝑔+𝑎 where [ · ]𝑗:𝑘 denotes the 𝑗th through 𝑘th rows of a vector or matrix. II.E.

Approximations

We shall now develop useful approximate expressions for the terms on the right hand sides of (10) and (11) as well as for the effects of (constant) control. The expressions will be based on a drag-only approximation for the acceleration contribution to 𝑎 in (2) caused by aerodynamic forces. This is motivated by the fact that, for a course corrected projectile in essentially ballistic flight (with 𝑢 = 0), the aerodynamic angles are small and the main contribution from the aerodynamic forces is (zero-lift) drag. It will be useful in a few places below to employ the notation 𝑥(𝑡) = [Φ((𝑥0 , 𝑣 0 ), 𝑡)]1:3 ,

𝑣(𝑡) = [Φ((𝑥0 , 𝑣 0 ), 𝑡)]4:6 ,

(13)

and accordingly set 𝑉 (𝑡) = ‖𝑣(𝑡)‖. II.E.1.

Subsonic drag-only approximation

The drag-only approximation for 𝑎 in (2) is defined by setting 𝑎((𝑥, 𝑣), 𝑡) = 𝑎𝑑 ((𝑥, 𝑣))

(14)

where

𝜌(𝑥)𝑆𝐶𝐷 𝑉 𝑣. (15) 2𝑚 Here, 𝜌 is the air density, 𝑆 is the reference area, 𝑚 is the mass, 𝐶𝐷 is the drag force aerodynamic coefficient and 𝑉 = ‖𝑣‖ is the air speed (neglecting wind). The right hand side in (15) is often expressed in terms of the ballistic coefficient 𝛽𝑑 = 𝑚/(𝑆𝐶𝐷 ), here defined in units of mass per area. The relation (15) gives 𝑎𝑑 ((𝑥, 𝑣)) = −

)︀𝑇 𝜕𝑎𝑑 ((𝑥, 𝑣)) 𝑆𝐶𝐷 𝑉 (︀ =− 𝑣 ∇𝜌(𝑥) . 𝜕𝑥 2𝑚 7 of 26 American Institute of Aeronautics and Astronautics Copyright © 2013 by John W.C. Robinson, Peter Strömbäck. Published by the American Institute of Aeronautics and Astronautics, Inc., with permission.

(16)

If we moreover focus on the subsonic part of the trajectory, which for almost all firing cases is synonymous with the part of the trajectory that follows after the initial deceleration, we may use the approximation that 𝐶𝐷 is constant.k With this approximation we have )︀ 𝜌(𝑥)𝑆𝐶𝐷 𝑉 (︀ 1 𝜕𝑎𝑑 ((𝑥, 𝑣)) =− 𝐼 3×3 + 2 𝑣𝑣 𝑇 . 𝜕𝑣 2𝑚 𝑉

(17)

Size estimates. A typical value of the air density at sea level is 1.225kg/m3 and the absolute value of the (vertical) air density gradient is less than 2.0×10−4 kg/m3 /m up to an altitude of 11000m (tropopause) in the ISA model. Hence, for a typical 155mm artillery shell with mass 𝑚 = 43kg and reference area 𝑆 = 0.189m2 traveling at an airspeed 𝑉 = 300m/s, and with a drag coefficient 𝐶𝐷 = 0.15 we have 𝑆𝐶𝐷 𝑉 = 0.01, 2𝑚

Downloaded by John Robinson on September 5, 2013 | http://arc.aiaa.org | DOI: 10.2514/6.2013-5114

and

(︀ )︀𝑇 0 ≤ |𝑣 ∇𝜌(𝑥) 𝜉| ≤ 0.06, )︀ (︀ 1 1 ≤ | 𝐼 3×3 + 2 𝑣𝑣 𝑇 𝜉| ≤ 2, 𝑉

(18)

𝜉 ∈ R3 , ‖𝜉‖ = 1, 𝜉 ∈ R3 , ‖𝜉‖ = 1.

(19) (20)

The ballistic coefficient 𝛽𝑑 is approximately 1.52 × 104 kg/m2 . II.E.2.

Effect of controls

¯ defined in (3). The solution Φ((𝑥0 , 𝑣 0 ), 𝑡; 𝑢) in (4) is nonlinear in the control 𝑢, which enters in the term 𝑎 However, the effect can often be well approximated as linear influence when 𝑢 is small in norm. Moreover, the expressions will be simple when 𝑢 can be finite dimensionally parametrized (such as polynomial, spline, etc.). The simplest case is when the control 𝑢 is constant in time, which is the case we will focus on. For a time invariant 𝑢 in we have, by 𝐶 1 -smoothness of 𝑎, Φ, thatl (︀ )︀ 𝜕𝑎 Φ((𝑥0 , 𝑣 0 ), 𝑡; 𝑢), 𝑡 ¯ ((𝑥0 , 𝑣 0 ), 𝑡; 𝑢)) 𝜕𝑎 = + 𝐼 3×3 = 𝜕𝑢 (︀ 𝜕𝑢 )︀ (︀ )︀ 𝜕𝑎 (𝑥, 𝑣(𝑡)), 𝑡 ⃒⃒ 𝜕𝑎 (𝑥(𝑡), 𝑣), 𝑡 ⃒⃒ 𝜕[Φ((𝑥0 , 𝑣 0 ), 𝑡; 𝑢)]1:3 𝜕[Φ((𝑥0 , 𝑣 0 ), 𝑡; 𝑢)]4:6 + + 𝐼 3×3 𝑥=𝑥(𝑡) 𝑣=𝑣(𝑡) 𝜕𝑥 𝜕𝑢 𝜕𝑣 𝜕𝑢

(21)

(with notation as in (12) and (13)). We shall be interested in the size and shape of the right hand side of (21) for time invariant 𝑢 when ‖𝑢‖ is small. In particular, for a vacuum trajectory (𝑎 = 0) we have ¯ ((𝑥0 , 𝑣 0 ), 𝑡; 0)) 𝜕𝑎 = 𝐼 3×3 . 𝜕𝑢

(22)

Drag-only approximation. To investigate more accurate approximations for the right hand side of (21) based on the drag-only approximation (15) we recall (4)–(6). From these relations it follows that, for time invariant 𝑢, we have ]︃ ∫︁ 𝑡 ∫︁ 𝑡 [︃ (︀ )︀ 𝜕 𝑎 ¯ ((𝑥0 , 𝑣 0 ), 𝑠; 0)) ¯ ((𝑥0 , 𝑣 0 ), 𝑠; 0)) 𝜕Φ((𝑥0 , 𝑣 0 ), 𝑡; 0) (𝑡 − 𝑠)𝐼 2×2 𝜕 𝑎 = exp (𝑡 − 𝑠)𝐴 𝐵 𝑑𝑠 = 𝑑𝑠. 𝜕𝑢 𝜕𝑢 𝜕𝑢 𝐼 2×2 0 0 (23) Combining this with (21) and (8) yields the functional equation (︀ )︀ ∫︁ 𝑡 𝜕𝑎 (𝑥, 𝑣(𝑡)), 𝑡 ⃒⃒ ¯ ((𝑥0 , 𝑣 0 ), 𝑡; 0)) ¯ ((𝑥0 , 𝑣 0 ), 𝑠; 0)) 𝜕𝑎 𝜕𝑎 = (𝑡 − 𝑠) 𝑑𝑠 𝑥=𝑥(𝑡) 𝜕𝑢 𝜕𝑥 𝜕𝑢 0 (︀ )︀ ∫︁ 𝑡 𝜕𝑎 (𝑥(𝑡), 𝑣), 𝑡 ⃒⃒ ¯ ((𝑥0 , 𝑣 0 ), 𝑠; 0)) 𝜕𝑎 + 𝑑𝑠 + 𝐼 3×3 . (24) 𝑣=𝑣(𝑡) 𝜕𝑣 𝜕𝑢 0 k This

is reasonable for subsonic flight, cf. e.g. Ref. 26, for small aerodynamic angles. of the solution Φ with respect to 𝑢 follows from standard properties of flows (continuity with respect to initial values) since 𝑢 here can be regarded as an extra state variable with zero time derivative. l Smoothness

8 of 26 American Institute of Aeronautics and Astronautics Copyright © 2013 by John W.C. Robinson, Peter Strömbäck. Published by the American Institute of Aeronautics and Astronautics, Inc., with permission.

¯ /𝜕𝑢 where the (integral) operator in question is contractive for This is an affine functional equationm for 𝜕 𝑎 sufficiently small 𝑡 and then the equation has a unique solution. With the drag-only approximation (15) and (16), (17) the equation (24) takes the formn ∫︁ (︀ )︀𝑇 𝑡 ¯ ((𝑥0 , 𝑣 0 ), 𝑠; 0)) ¯ ((𝑥0 , 𝑣 0 ), 𝑡; 0)) 𝑆𝐶𝐷 𝑉 (𝑡) 𝜕𝑎 𝜕𝑎 = 𝐼 3×3 − 𝑣(𝑡) ∇𝜌(𝑥(𝑡)) (𝑡 − 𝑠) 𝑑𝑠 𝜕𝑢 2𝑚 𝜕𝑢 0 ∫︁ )︀ 𝑡 𝜕 𝑎 ¯ ((𝑥0 , 𝑣 0 ), 𝑠; 0)) 1 𝜌(𝑥(𝑡))𝑆𝐶𝐷 𝑉 (𝑡) (︀ 𝑇 𝐼 3×3 + 𝑣(𝑡)𝑣(𝑡) − 𝑑𝑠 (25) 2𝑚 𝑉 (𝑡)2 𝜕𝑢 0 and we note that the vacuum trajectory solution (22) results if we let the ballistic coefficient tend to infinity. If we use the the vacuum trajectory solution (22) as a zeroth order approximation for the integrands on the right in (25) we obtain the first order approximation

Downloaded by John Robinson on September 5, 2013 | http://arc.aiaa.org | DOI: 10.2514/6.2013-5114

)︀ 𝑡2 𝑆𝐶𝐷 𝑉 (𝑡) (︀ )︀𝑇 ¯ ((𝑥0 , 𝑣 0 ), 𝑡; 0)) 𝜕𝑎 𝜌(𝑥(𝑡))𝑆𝐶𝐷 𝑉 (𝑡) (︀ 1 𝑇 𝑣(𝑡)𝑣(𝑡) = 𝐼 3×3 −𝑡 𝐼 3×3 + − 𝑣(𝑡) ∇𝜌(𝑥(𝑡)) . (26) 𝜕𝑢 2𝑚 𝑉 (𝑡)2 2 2𝑚 The size estimates in (18)–(20) indicate that for times 𝑡 up to about 10s the nonconstant terms on the right in (26) can be neglected, with an error of at most a few percent. Based on this we see that it can in fact be reasonableo more generally to use (22) as an approximation (i.e. use the zeroth order approximation), at least for short prediction times. Better approximations can be obtained from higher order approximations resulting from further iterations with (25). If we use (22) as an approximation we obtain from (23) that ]︃ [︃ 2 𝑡 𝜕Φ((𝑥0 , 𝑣 0 ), 𝑡; 0) 𝐼 3×3 2 , = 𝜕𝑢 𝑡𝐼 3×3 and a first order approximation of the effect of a time invariant 𝑢 in Φ in (4) can thus be obtained as Φ((𝑥0 , 𝑣 0 ), 𝑡; 𝑢) = Φ((𝑥0 , 𝑣 0 ), 𝑡; 0) +

𝜕Φ((𝑥0 , 𝑣 0 ), 𝑡; 0) 𝑢 + 𝑜(‖𝑢‖) = 𝜕𝑢 Φ((𝑥0 , 𝑣 0 ), 𝑡) +

(︀ 𝑡2 )︀ 𝑢, 𝑡𝑢 + 𝑜(‖𝑢‖). (27) 2

This is the basic relation for the influence of control that we shall use. II.E.3.

Acceleration evolution

In Appendix VI several forms of approximate expressions for the right hand sides of (10) and (11) are described based on the drag-only approximation (15) and (16), (17). In particular, it is shown that the upper part of the matrix on the right in (10) can be approximated as ∫︁

𝑡

(𝑡 − 𝑠)

𝐼 3×3 + 0

)︀𝑇 ¯ ((𝑥, 𝑣 0 ), 𝑠) ⃒⃒ 𝜕𝑎 𝑡2 𝑆𝐶𝐷 𝑉0 (︀ 𝑣 𝑑𝑠 = 𝐼 0 ∇𝜌(𝑥0 ) 3×3 − 𝑥=𝑥 0 𝜕𝑥 2 2𝑚

(28)

and for the lower part we have the approximation ∫︁ 𝑡 )︀ (︀ )︀𝑇 )︀𝑇 𝑡2 (︀ 𝑆𝐶𝐷 𝑉0 )︀2 (︀ ¯ ((𝑥, 𝑣 0 ), 𝑠) ⃒⃒ 𝜕𝑎 𝑆𝐶𝐷 𝑉0 (︀ 1 𝑑𝑠 = −𝑡 𝑣 0 ∇𝜌(𝑥0 ) + 𝜌(𝑥0 ) 𝐼 3×3 + 2 𝑣 0 𝑣 𝑇0 𝑣 0 ∇𝜌(𝑥0 ) . 𝑥=𝑥 0 𝜕𝑥 2𝑚 2 2𝑚 𝑉0 0 (29) For the upper part of the matrix on the right in (11) we have the approximation ∫︁

𝑡

(𝑡 − 𝑠)

𝑡𝐼 3×3 + 0

)︀ ¯ ((𝑥0 , 𝑣), 𝑠) ⃒⃒ 𝜕𝑎 𝑡2 𝜌(𝑥0 )𝑆𝐶𝐷 𝑉0 (︀ 1 𝑑𝑠 = 𝑡𝐼 − 𝐼 3×3 + 2 𝑣 0 𝑣 𝑇0 3×3 𝑣=𝑣 0 𝜕𝑣 2 2𝑚 𝑉0

m It

(30)

can be interpreted as an inhomogeneous Volterra equation of the second kind. this form, it is easy to see that the functional equation has unique solution when the ballistic coeffient is high enough. o Neglecting the higher order terms is equivalent to introducing an “matched disturbance” in the parlance of robust control.29

n In

9 of 26 American Institute of Aeronautics and Astronautics Copyright © 2013 by John W.C. Robinson, Peter Strömbäck. Published by the American Institute of Aeronautics and Astronautics, Inc., with permission.

and for the lower part we have ∫︁ 𝐼 3×3 +

Downloaded by John Robinson on September 5, 2013 | http://arc.aiaa.org | DOI: 10.2514/6.2013-5114

0

𝑡

)︀ ¯ ((𝑥0 , 𝑣), 𝑠) ⃒⃒ 𝜕𝑎 𝜌(𝑥0 )𝑆𝐶𝐷 𝑉0 (︀ 1 𝑑𝑠 = 𝐼 3×3 − 𝑡 𝐼 3×3 + 2 𝑣 0 𝑣 𝑇0 𝑣=𝑣 0 𝜕𝑣 2𝑚 𝑉0 (︁ 2 )︀2 )︁ (︀ )︀ 𝑆𝐶𝐷 𝑉0 (︀ 1 𝑡 𝑆𝐶𝐷 𝑉0 𝑇 . (31) 𝑣 0 ∇𝜌(𝑥0 ) − 𝜌(𝑥0 )2 𝐼 3×3 + 2 𝑣 0 𝑣 𝑇0 − 2 2𝑚 2𝑚 𝑉0

It is clear from all these expressions that as the ballistic coefficient tends to infinity only the vacuum trajectory terms remain and that these terms generally provide a good approximation for small 𝑡, as alluded to above. To see when the various terms come into play, recall the size estimates (18)–(20). Applied to (28), these estimates suggest that for times 𝑡 larger that about 10s the term that is quadratic in 𝑡 needs to be taken into account, if an acceptable relative error is in the order of a few percent. For (29) the size estimates similarly suggest that for times 𝑡 larger that about a few seconds the term that is quadratic in 𝑡 needs to be taken into account, for a relative error of a few percent. Similarly, for (30) the size estimates (18)–(20) suggest that for times 𝑡 that are larger than about a couple of seconds the term that is quadratic in 𝑡 need to be taken into account, for a relative error of a few percent. Finally, for (31) the size estimates suggest that for times larger than about one second the term that is linear in 𝑡 needs to be taken into account, if an acceptable relative error is in the order of a few percent. The term that is quadratic in 𝑡 needs to be taken into account for times larger than about 5s. However, we have here only considered how many terms need to be taken into account for a certain relative accuracy with respect to the approximation itself, not with respect to the underlying true quantity. Moreover, equal relative accuracy is not always needed when approximating the integrals in (10) and (11); it all depends on the typical sizes (and directions) of the perturbations ∆𝑥 and ∆𝑣 in (9) and on the prediction horizon.p

III.

Guidance

We shall now complement the perturbation expressions (9)–(11) (and associated approximations in (28)– (31)) and the control influence relation (27) with explicit equations for guidance. More precisely, we shall derive explicit expressions for the acceleration commands to be issued by guidance to the autopilot function of the course correcting fuze. III.A.

Predictive guidance

Consider a firing solution flow Φ𝑥𝐺 with a firing solution (𝑥𝐺 , 𝑥𝑇 , 𝑣 𝐺 ) and nominal trajectory 𝑡 ↦→ Φ𝑥𝐺 ((𝑥𝐺 , 𝑣 𝐺 ), 𝑡),

(32)

Φ𝑥𝐺 ((𝑥𝐺 , 𝑣 𝐺 ), 𝑡𝑇 ) = (𝑥𝑇 , 𝑣 𝑇 ).

(33)

with terminal condition For such a firing solution we have the (flow) group property (︀ )︀ Φ𝑥𝐺 Φ𝑥𝐺 ((𝑥𝐺 , 𝑣 𝐺 ), 𝑡1 ), 𝑡𝑇 − 𝑡1 = (𝑥𝑇 , 𝑣 𝑇 ),

(34)

where 𝑡1 ∈ [0, 𝑡𝑇 ] is any time point in the nominal time-of-flight interval. Using this property, the deviation in terminal state obtained when starting at time 𝑡1 at some state (𝑥1 , 𝑣 1 ), near the nominal trajectory, can conveniently be calculated (approximately) using the perturbation approximations developed above. III.A.1.

Zero effort terminal deviation

Let 𝑡1 ∈ [0, 𝑡𝑇 ] ⊂ ℐ and (𝑥1 , 𝑣 1 ) ∈ ℳ be an arbitrary initial time and initial state, respectively, for a trajectory in the firing solution flow which is a perturbation of the nominal trajectory. Then, using the terminal condition (33) and the group property (34), the deviation 𝜁 0 (𝑡𝐹 ) ∈ R6 between the terminal state p Note also that in the guidance algorithm in section III.A.3 below only (28) and (30) (or their equivalents, obtained by other approximation methods) are used.

10 of 26 American Institute of Aeronautics and Astronautics Copyright © 2013 by John W.C. Robinson, Peter Strömbäck. Published by the American Institute of Aeronautics and Astronautics, Inc., with permission.

Φ𝑥𝐺 ((𝑥1 , 𝑣 1 ), 𝑡𝐹 − 𝑡1 ) at time 𝑡𝐹 ∈ ℐ of the perturbed trajectory and the nominal terminal state (𝑥𝑇 , 𝑣 𝑇 ) can be written 𝜁 0 (𝑡𝐹 ) = Φ𝑥𝐺 ((𝑥1 , 𝑣 1 ), 𝑡𝐹 − 𝑡1 ) − (𝑥𝑇 , 𝑣 𝑇 ) = Φ𝑥𝐺 ((𝑥1 , 𝑣 1 ), 𝑡𝐹 − 𝑡1 ) − Φ𝑥𝐺 (Φ𝑥𝐺 ((𝑥𝐺 , 𝑣 𝐺 ), 𝑡1 ), 𝑡𝑇 − 𝑡1 ). (35) The deviation 𝜁 0 (𝑡𝐹 ) in (35) will be called the zero effort (terminal) deviation (at time 𝑡𝐹 ). For the ZED we can apply the first order perturbation approximation in (9) to obtain, approximately (for small ∆𝑥, ∆𝑣, ∆𝑡) 𝜕Φ𝑥𝐺 (Φ𝑥𝐺 ((𝑥𝐺 , 𝑣 𝐺 ), 𝑡1 ), 𝑡𝑇 − 𝑡1 )∆𝑥 𝜕𝑥 𝜕Φ𝑥𝐺 𝜕Φ𝑥𝐺 (Φ𝑥𝐺 ((𝑥𝐺 , 𝑣 𝐺 ), 𝑡1 ), 𝑡𝑇 − 𝑡1 )∆𝑣 + (Φ𝑥𝐺 ((𝑥𝐺 , 𝑣 𝐺 ), 𝑡1 ), 𝑡𝑇 − 𝑡1 )∆𝑡, + 𝜕𝑣 𝜕𝑡

𝜁 0 (𝑡𝐹 ) = Φ𝑥𝐺 ((𝑥1 , 𝑣 1 ), 𝑡𝐹 − 𝑡1 ) − (𝑥𝑇 , 𝑣 𝑇 ) =

(36)

where ∆𝑥 = 𝑥1 − [Φ𝑥𝐺 ((𝑥𝐺 , 𝑣 𝐺 ), 𝑡1 )]1:3 ,

Downloaded by John Robinson on September 5, 2013 | http://arc.aiaa.org | DOI: 10.2514/6.2013-5114

III.A.2.

∆𝑣 = 𝑣 1 − [Φ𝑥𝐺 ((𝑥𝐺 , 𝑣 𝐺 ), 𝑡1 )]4:6 ,

∆𝑡 = 𝑡𝐹 − 𝑡𝑇 .

Terminal position deviation with control

It is straightforward to compute an approximation to the deviation in terminal state obtained when a (small) constant control acceleration 𝑢 is applied to the firing flow using the relation (27) and the ZED expressions above. With notation as in section III.A.1, consider the deviation 𝜁 𝑢 (𝑡𝐹 ) = Φ𝑥𝐺 ((𝑥1 , 𝑣 1 ), 𝑡𝐹 − 𝑡1 ; 𝑢) − (𝑥𝑇 , 𝑣 𝑇 ) obtained when a constant control acceleration 𝑢 is applied to the perturbed trajectory. From (27) we then have that the deviation in terminal state with constant control is approximately (for small 𝑢) 𝜁 𝑢 (𝑡𝐹 ) = Φ𝑥𝐺 ((𝑥1 , 𝑣 1 ), 𝑡𝐹 − 𝑡1 ; 𝑢) − (𝑥𝑇 , 𝑣 𝑇 ) = )︀ (︀ (𝑡𝐹 − 𝑡1 )2 )︀ (︀ (𝑡𝐹 − 𝑡1 )2 𝑢, (𝑡𝐹 − 𝑡1 )𝑢 = 𝜁 0 (𝑡𝐹 ) + 𝑢, (𝑡𝐹 − 𝑡1 )𝑢 , (37) Φ𝑥𝐺 ((𝑥1 , 𝑣 1 ), 𝑡𝐹 − 𝑡1 ) − (𝑥𝑇 , 𝑣 𝑇 ) + 2 2 where 𝜁 0 (𝑡𝐹 ) is the ZED in (35). Since the state vector (𝑥, 𝑣) is six dimensional but the control 𝑢 is only three dimensional there is in general no solution for the problem of driving the entire state deviation in (37) to zero during [𝑡𝐹 − 𝑡1 , 𝑡𝑇 ] using a constant control acceleration 𝑢 (regardless of 𝑡𝐹 ). Thus, some form of trade-off must be made, e.g. by only considering a condition on the position part of the state, or by introducing a weighted form of miss distance to be minimized. In the course corrected projectile application the former alternative is the relevant one and we have approximately (for small 𝑢), using (37), that [𝜁 𝑢 (𝑡𝐹 )]1:3 = [Φ𝑥𝐺 ((𝑥1 , 𝑣 1 ), 𝑡𝐹 − 𝑡1 ; 𝑢)]1:3 − 𝑥𝑇 = [Φ𝑥𝐺 ((𝑥1 , 𝑣 1 ), 𝑡𝐹 − 𝑡1 )]1:3 − 𝑥𝑇 +

(𝑡𝐹 − 𝑡1 )2 (𝑡𝐹 − 𝑡1 )2 𝑢 = 𝜁 0,𝑥 (𝑡𝐹 ) + 𝑢 2 2

(38)

where 𝜁 0,𝑥 (𝑡𝐹 ) ∈ R3 is the zero effort (terminal) position deviation (at time 𝑡𝐹 ) (ZEPD) defined by 𝜁 0,𝑥 (𝑡𝐹 ) = [𝜁 0 (𝑡𝐹 )]1:3 = [Φ𝑥𝐺 ((𝑥1 , 𝑣 1 ), 𝑡𝐹 − 𝑡1 )]1:3 − 𝑥𝑇 .

(39)

Using the ZEPD, the (relative) closest point of approach (CPA) 𝜁 0,𝑥 (𝑡𝐶𝑃 𝐴 ) ∈ R3 on a perturbed trajectory, relative to the nominal terminal position 𝑥𝑇 , can be determined. The CPA is defined by the condition that for some 𝑡𝐶𝑃 𝐴 ∈ ℐ it holds that ‖𝜁 0,𝑥 (𝑡𝐹 )‖ ≥ ‖𝜁 0,𝑥 (𝑡𝐶𝑃 𝐴 )‖,

𝑡𝐹 ∈ ℐ.

(40)

Determining the CPA is an instance of the point projection problem for a curve.30 A necessary condition for the inequality (40) to hold is that the the orthogonality condition (“footpoint condition”) 𝜁 0,𝑥 (𝑡𝐶𝑃 𝐴 )𝑇 [Φ𝑥𝐺 ((𝑥1 , 𝑣 1 ), 𝑡𝐶𝑃 𝐴 − 𝑡1 )]4:6 = 0

(41)

is satisfied. If the perturbation of the nominal trajectory is small enough there will be a unique 𝑡𝐶𝑃 𝐴 such that this holds and we shall assume that this is the case. The quantity 𝑑𝑍𝐸𝑀 = ‖𝜁 0,𝑥 (𝑡𝐶𝑃 𝐴 )‖ will be called the zero effort miss (distance) (ZEM).

11 of 26 American Institute of Aeronautics and Astronautics Copyright © 2013 by John W.C. Robinson, Peter Strömbäck. Published by the American Institute of Aeronautics and Astronautics, Inc., with permission.

By combining the perturbation expression (9) and the derivative (12)) we can obtain an approximation for the position [Φ𝑥𝐺 ((𝑥1 , 𝑣 1 ), 𝑡𝐹 − 𝑡1 )]1:3 (at an arbitrary terminal time 𝑡𝐹 ) as [Φ𝑥𝐺 ((𝑥1 , 𝑣 1 ), 𝑡𝐹 − 𝑡1 )]1:3 = [Φ𝑥𝐺 ((𝑥1 , 𝑣 1 ), 𝑡𝑇 − 𝑡1 )]1:3 + (𝑡𝐹 − 𝑡𝑇 )[Φ𝑥𝐺 ((𝑥1 , 𝑣 1 ), 𝑡𝑇 − 𝑡1 )]4:6 , and an associated estimate for 𝜁 0,𝑥 (𝑡𝐶𝑃 𝐴 ) follows from (39) as 𝜁 0,𝑥 (𝑡𝐶𝑃 𝐴 ) = 𝜁 0,𝑥 (𝑡𝑇 ) + (𝑡𝐶𝑃 𝐴 − 𝑡𝑇 )[Φ𝑥𝐺 ((𝑥1 , 𝑣 1 ), 𝑡𝑇 − 𝑡1 )]4:6 .

(42)

If we further use the approximations [Φ𝑥𝐺 ((𝑥1 , 𝑣 1 ), 𝑡𝐶𝑃 𝐴 − 𝑡1 )]4:6 = [Φ𝑥𝐺 ((𝑥1 , 𝑣 1 ), 𝑡𝑇 − 𝑡1 )]4:6 = 𝑣 𝑇 , we can solve for 𝑡𝐶𝑃 𝐴 − 𝑡𝑇 and 𝜁 0,𝑥 (𝑡𝐶𝑃 𝐴 ) using (41), (42). This yields the approximate relations

Downloaded by John Robinson on September 5, 2013 | http://arc.aiaa.org | DOI: 10.2514/6.2013-5114

𝑡𝐶𝑃 𝐴 − 𝑡𝑇 = −

𝜁 0,𝑥 (𝑡𝑇 )𝑇 𝑣 𝑇 ‖𝑣 𝑇 ‖2

(︀ 𝑣 𝑇 𝑣 𝑇𝑇 )︀ and 𝜁 0,𝑥 (𝑡𝐶𝑃 𝐴 ) = 𝐼 3×3 − 𝜁 (𝑡𝑇 ) = 𝑃 [𝑣𝑇 ]⊥ 𝜁 0,𝑥 (𝑡𝑇 ), ‖𝑣 𝑇 ‖2 0,𝑥

(43)

where 𝑃 [𝑣𝑇 ]⊥ denotes projection onto the subspace orthogonal to 𝑣 𝑇 . Driving the ZEM to zero. If we take 𝑡𝐹 = 𝑡𝐶𝑃 𝐴 in (38) we obtain a simple approximate relation (for small 𝑢) that can be used to chose a constant control acceleration 𝑢 to drive the predicted terminal position ˆ where deviation to (approximately) zero during [𝑡1 , 𝑡𝐶𝑃 𝐴 − 𝑡1 ]. The solution is to set 𝑢 = 𝑢 ˆ =− 𝑢

2 𝜁 (𝑡𝐶𝑃 𝐴 ), (𝑡𝐶𝑃 𝐴 − 𝑡1 )2 0,𝑥

(44)

and 𝜁 0,𝑥 (𝑡𝐶𝑃 𝐴 ) is the (uncontrolled) CPA in (40), which can be approximated as in (43). Useful approximations to 𝜁 0,𝑥 (𝑡𝑇 ) (the ZEPD at the nominal terminal time 𝑡𝑇 ) occurring on the right in (43) can be obtained from (9)–(11) and the expressions in sections II.E.3 and III.A.2 (cf. also section VI.D in appendix). III.A.3.

Guidance algorithm

The resulting predictive guidance algorithm can be described as follows. Before flight, compute a nominal trajectory (32) using any kind of high fidelity method. Then, at each guidance update instant during flight, do the following: (1) Compute an approximation to the ZEPD 𝜁 0,𝑥 (𝑡𝐹 ) in (39) for 𝑡𝐹 = 𝑡𝑇 (the nominal terminal time) using the current position 𝑥1 , velocity 𝑣 1 , and time 𝑡1 . The quantity [Φ𝑥𝐺 ((𝑥1 , 𝑣 1 ), 𝑡𝑇 − 𝑡1 )]1:3 is then approximated by setting 𝑥1 = 𝑥0 + ∆𝑥, 𝑣 1 = 𝑣 0 + ∆𝑣, 𝑡 = 𝑡𝑇 − 𝑡1 and considering the first three terms on the right in (9). The Jacobians on the right in (9) are obtained from (10), (11) where the integrals can be approximated as in section II.E.3 (or as in Appendix). (2) Compute approximations of 𝑡𝐶𝑃 𝐴 − 𝑡𝑇 and 𝜁 0,𝑥 (𝑡𝐶𝑃 𝐴 ) in (40) using the formulas in (43). (3) Compute the acceleration command using (44). The guidance commands must then be translated into control actuator settings, for instance using a guidance command filter and autopilot as described in Appendix. III.B.

Path Following Guidance

A path following guidance principle can be obtained from the predictive guidance principle devised above by simply reinterpreting the terminal time 𝑡𝑇 and state (𝑥𝑇 , 𝑣 𝑇 ) as a state a bit further down on the nominal trajectory (rather than the terminal state).

IV.

Simulations

The theory derived in the previous sections shall now be illustrated with simulations utilizing a full 6DOF nonlinear model of a spinning 155mm spinning artillery shell with a generic 2D course correcting fuze. 12 of 26 American Institute of Aeronautics and Astronautics Copyright © 2013 by John W.C. Robinson, Peter Strömbäck. Published by the American Institute of Aeronautics and Astronautics, Inc., with permission.

IV.A.

Simulation model

The simulations are performed using a full nonlinear 6DOF model with equations of motions formulated in the non-spinning frame31 with a flat Earth gravity model and a standard ICAO atmosphere model. A 4th order Runge-Kutta integration method with fixed time step of 0.001s is used to solve the equations of motion. Effects of Earth rotation are neglected in the simulations. The guidance commands are computed using the steps outlined in section III.A.3 above and the guidance loop operates in continuous time. An autopilot, described in appendix, calculates trim values for the aerodynamic angles and angular rates which are used to translate the guidance commands to control surface deflection forces. The autopilot does not include any dynamics syntheses but instead a simple filter, also this described in appendix, is inserted to act on the guidance commands before these are passed to the autopilot. IV.A.1.

Shell properties

Downloaded by John Robinson on September 5, 2013 | http://arc.aiaa.org | DOI: 10.2514/6.2013-5114

The simulated shell is a model of a M107 artillery shell with mass properties, dimensions and aerodynamic coefficients obtained from Ref. 32. The unperturbed initial velocity of the projectile as it leaves the muzzle is 684.3m/s and the initial spin rate is 1102.6rad/s. The shell data is summarized in table 1. Table 1: Properties of the simulated shell. Description Total mass Caliber Axial moment of inertia Lateral moment of inertia Muzzle velocity Initial spin rate

IV.A.2.

Value 43.0kg 0.155m 0.144kg m2 1.216kg m2 684.3m/s 1102.6rad/s (175.45 rps)

Fuze guidance and control properties

The fuze is modeled as a generic 2D course correction fuze where control forces are applied on the nose of the projectile, representing the effect of normal force from deflectable canards. The distance from the center of gravity to the center of pressure for the canards was set to 0.4m. The data for the available control force from canards was taken from Ref. 12. A limit of 80N for the canard control force was used for all simulations except those relating to the influence of the control force limit. The value 80N is slightly over the value considered in Ref. 12 but is well within the achievable range with respect to packaging and stability requirements.12 Control force limits and guidance filter parameters are summarized in table 2 (the notation is explained in section VII in Appendix where the guidance command filter is described in more detail). In the guidance calculations, the integrals on the right in (10), (11) (cf. step Table 2: Properties and limitations of the guidance and control model. Description Commanded acceleration limit (𝑐𝑙 ) Commanded acceleration rate limit (𝑟𝑙 ) Canard control force limit Guidance filter natural frequency (𝜔𝑔 ) Guidance filter damping (𝜁𝑔 ) Canard control force moment arm

Value ±20m/s2 ±0.02 · 9.81m/s2 /s ±80N 2 · 2𝜋rad/s 0.8 0.4m

1 in the guidance algorithm in section III.A.3) were approximated as in (28) and (30). The approximations obtained from (28) and (30) when retaining only terms only constant and linear in 𝑡 will be called the first order approximation. The first order approximation provides the same corrections as for a vacuum trajectory.

13 of 26 American Institute of Aeronautics and Astronautics Copyright © 2013 by John W.C. Robinson, Peter Strömbäck. Published by the American Institute of Aeronautics and Astronautics, Inc., with permission.

IV.B.

Results

The results from simulations of the guided and unguided artillery shell will be presented below where the chosen figure of merit is mean value and standard deviation of miss distance, downrange and cross range relative to a nominal unperturbed trajectory. The miss distance is defined as the distance 𝑑 = ‖x0 − x𝑘 ‖, between the actual ground impact point x𝑘 and the nominal (unperturbed) impact point x0 . IV.B.1.

Nominal trajectory

Downloaded by John Robinson on September 5, 2013 | http://arc.aiaa.org | DOI: 10.2514/6.2013-5114

A nominal firing case was chosen around which perturbations were studied. For the nominal unperturbed firing case an initial gun elevation of 30.0∘ and muzzle velocity of 684.3m/s was chosen. The trajectory, total angle of attack and velocity profile are presented in figures 1 and 2. The simulation of the (unguided) shell results in a ground impact after 49.08s at 15890.89m downrange and 296.35m crossrange with apogee at 3160.39m. During most of the flight the shell travels at subsonic speeds. The spin rate is damped from the initial 1102.6rad/s to 829.6rad/s at impact.

(a)

(b)

Figure 1: Trajectory for a M107 projectile launched at an initial elevation of 30∘ .

(a) Total angle of attack.

(b) Velocity.

Figure 2: Total angle of attack and velocity of the projectile over time. The dashed line represents passing Mach 1.

14 of 26 American Institute of Aeronautics and Astronautics Copyright © 2013 by John W.C. Robinson, Peter Strömbäck. Published by the American Institute of Aeronautics and Astronautics, Inc., with permission.

Downloaded by John Robinson on September 5, 2013 | http://arc.aiaa.org | DOI: 10.2514/6.2013-5114

A plot of the root locus for the linearized dynamics around free flight trim and the associated damping and natural (angular) frequency the two pole pairs are given in figure 3. The nutation (fast) mode is weakly damped throughout the flight, with a minimum near apogee, and the precession (slow) mode becomes poorly damped on the downleg. Since the natural frequency of the precession is in the order of 1Hz it has the potential of affecting the dynamics in the closed loop formed when guidance is engaged.

(a) Root locus plot.

(b) Damping (𝜁) and natural frequency (𝜔0 ).

Figure 3: Root locus (left) and damping 𝜁 and natural frequency 𝜔0 (right) of the precession (green) and nutation (blue) modes. The minimal values for the damping of the two modes are indicated.

IV.B.2.

Dispersion analysis

The capability to handle initial errors are studied by comparing unguided and guided dispersions for individual errors and combinations of errors. The nominal firing case described above is used as reference and initial errors in elevation, azimuth and velocity are added as perturbations. The elevation and azimuth errors are normally distributed with a standard deviation of 0.1∘ . The muzzle velocity errorq is also normally distributed but with a standard deviation of 3.0m/s. For each studied set of error combinations, 1000 Monte Carlo simulations are performed and presented below as dispersion plots and tables describing the statistics. The results for standalone elevation and velocity errors are presented in figure 4a and 4b with statistics presented in table 3. The combinations of errors are presented in figure 5a and 5b with statistics in table 4. One of the largest errors affecting the dispersion of artillery shells is the variation in muzzle velocity. The results from perturbations in initial muzzle velocity presented in figure 4b and table 3 show a decrease in mean of miss distance from almost 59m unguided to a mean miss distance of slightly above 5m with the proposed 2D guidance. The corresponding miss distance standard deviation decreases from about 45m to about 10m. Note also that for all examined cases of perturbations the mean downrange error is slightly biased for the guided shell which travels a few meters shorter than the unguided shell, on average. From table 4 it can be seen that introducing azimuth errors (as compared to only combination of elevation and velocity errors) degrades the guided solution only slightly. The standard deviation of downrange errors is approximately 78m unguided for both cases and is reduced to close to 12m with guidance. The crossrange error standard deviation increases from approximately 3m to 28m when introducing azimuth errors to the unguided shell. The guided shell is able to decrease the dispersion to a standard deviation of close to 1m. In figure 6a and 6b the miss distance distribution corresponding to the dispersion plot in figure 5b is presented. A resemblance with the Rayleigh distribution can be noted. In the dispersion simulations a canard control force limit of 80N was used. The effect of saturated controls can be seen in figure 4b where the guided fuze is not always able to adequately reduce the effect q Both the mean and standard deviation of the muzzle velocity error vary considerably with gun wear33 and we have here chosen a figure for the standard deviation based on the results in Ref. 34.

15 of 26 American Institute of Aeronautics and Astronautics Copyright © 2013 by John W.C. Robinson, Peter Strömbäck. Published by the American Institute of Aeronautics and Astronautics, Inc., with permission.

Table 3: Statistics from 1000 Monte Carlo simulations where perturbations on the initial elevation and velocity have been added.

Miss distance [m] Downrange [m] Crossrange [m]

Elevation perturbation Unguided Guided Mean Std. Dev. Mean Std. Dev. 17.726 13.034 0.854 0.504 0.236 21.955 -0.348 0.862 0.022 1.530 0.061 0.341

Velocity perturbation Unguided Guided Mean Std. Dev. Mean Std. Dev. 58.755 44.621 5.393 9.521 3.358 73.689 -2.704 10.586 0.110 2.303 -0.112 0.622

Downloaded by John Robinson on September 5, 2013 | http://arc.aiaa.org | DOI: 10.2514/6.2013-5114

Table 4: Statistics from 1000 Monte Carlo simulations where perturbations on the initial elevation, velocity and azimuth have been added.

Miss distance [m] Downrange [m] Crossrange [m]

Elevation and velocity perturbation Unguided Guided Mean Std. Dev. Mean Std. Dev. 61.917 48.041 5.995 10.802 3.599 78.259 -3.180 11.917 0.133 2.848 -0.122 0.706

Elevation, velocity Unguided Mean Std. Dev. 70.078 44.950 3.552 78.261 1.288 28.236

and azimuth pert. Guided Mean Std. Dev. 6.120 11.021 -3.323 12.123 -0.119 0.956

of the initial velocity error. Since the added velocity errors are normally distributed it is fully possible that some of the errors are significantly higher than the standard deviation of 3.0m. The largest velocity errors observed in data for figure 4b was 10.3m/s. In figure 7a and 7b the commanded acceleration and applied control forces as function of time are visualized for this extreme case. Even though such a large error leads to saturated control in the vertical direction for a considerable amount of the flight time the error is still reduced significantly compared to the unguided shell. It should also be noted that even though the proposed 2D guidance performs well in terms of dispersion the orientation of the projectile at impact will be affected. This effect is not studied in this paper. In all presented cases of perturbations of initial errors the proposed predictive guidance significantly reduced the dispersion in downrange and crossrange for the studied projectile. It should be noted that the simulated shell is stable throughout the nominal trajectory as can be seen by the root locus plot of the linearized dynamics in figure 3. Shells with other aerodynamics and stability properties will have different dispersion results. In this section only a small subset of initial errors were studied. In reality variations in projectile properties such as aerodynamics coefficients, wind gusts and mass distribution occur. IV.B.3.

Miss distance as a function of control force limitations

In the previous section the dispersion for a guided shell with a control force limited to 80N was studied. In figure 8a and 8b miss distance is displayed as a function of the limit in control force. The maximal allowed control forcer is varied from 10N up to 140N for the same scenario as above with perturbations in initial velocity, elevation and azimuth with standard deviations of 3.0m/s, 0.1∘ and 0.1∘ respectively. This corresponds to the case presented in figure 5b and table 4. In figure 8b it can be seen that even for smaller control forces the miss distance (both mean value and standard deviation) is reduced significantly compared to the unguided shell. For a limit of 30N the mean miss distance is approximately 20m which is still much less than the unguided 70m.

r The

figure 140N is the maximal possible for the 155mm canard controlled shell studied in Ref. 12 before stability is lost.

16 of 26 American Institute of Aeronautics and Astronautics Copyright © 2013 by John W.C. Robinson, Peter Strömbäck. Published by the American Institute of Aeronautics and Astronautics, Inc., with permission.

(a) Perturbations in elevation only.

(b) Perturbations in velocity only.

Downloaded by John Robinson on September 5, 2013 | http://arc.aiaa.org | DOI: 10.2514/6.2013-5114

Figure 4: Dispersions for a guided and unguided projectile. Black ellipses represent 1𝜎 standard deviation.

(a) Perturbations in velocity and elevation.

(b) Perturbations in velocity, elevation and azimuth.

Figure 5: Dispersions for a guided and unguided projectile. Black ellipses represent 1𝜎 standard deviation.

(a) Miss distance probability density.

(b) Miss distance cum. probability distribution.

Figure 6: Miss distance probability density and cumulative distribution for the perturbation of initial values of elevation, velocity and azimuth.

17 of 26 American Institute of Aeronautics and Astronautics Copyright © 2013 by John W.C. Robinson, Peter Strömbäck. Published by the American Institute of Aeronautics and Astronautics, Inc., with permission.

Downloaded by John Robinson on September 5, 2013 | http://arc.aiaa.org | DOI: 10.2514/6.2013-5114

(a) Commanded acceleration in the fire control frame. (Note that saturation of commanded acceleration is performed in the non-spinning body frame.)

(b) Applied control forces in the non-spinning body frame.

Figure 7: Commanded acceleration and saturated control forces for the case of an initial velocity error of 10.3m/s.

(a) Mean of miss distance.

(b) Standard deviation of miss distance.

Figure 8: Miss distance as a function of maximal canard control force. Each data point represents the mean value or standard deviation from 100 Monte Carlo simulations.

18 of 26 American Institute of Aeronautics and Astronautics Copyright © 2013 by John W.C. Robinson, Peter Strömbäck. Published by the American Institute of Aeronautics and Astronautics, Inc., with permission.

IV.B.4.

Miss distance as a function of elevation

Downloaded by John Robinson on September 5, 2013 | http://arc.aiaa.org | DOI: 10.2514/6.2013-5114

To further analyze the expected accuracy of the proposed predictive control method the miss distance of the projectile was evaluated for a set of nominal elevations. In figure 9a the mean value of the miss distance is presented based on 100 Monte Carlo runs for each elevation. The initial velocity, elevation and azimuth errors are normally distributed with a standard deviation of 3.0m/s, 0.1∘ and 0.1∘ respectively. For very low elevations the effectiveness of the guidance is degraded. This can also be seen in figure 10 where the impact of the individual errors on the accuracy are displayed. For low elevation firing cases errors in the initial elevation is the most severe for the guided shell. For elevations above 10∘ the largest remaining error source is the initial velocity error. The elevation and azimuth error is however almost completely eliminated.

(a) Guided and an unguided projectile.

(b) First and higher order guidance terms (sec. IV.A.2).

Figure 9: Mean value of miss distance as a function of elevation.

(a) Unguided.

(b) Guided.

Figure 10: Miss distance standard deviation as a function of elevation for independent errors. In figure 9b the effect of using only first order terms in the guidance calculation (cf. sec. IV.A.2) is shown for a same set of elevations as shown in figure 9a. The results indicate that the here used estimates of the higher order terms only have a small impact on the miss distance and that they are not always accurate. In fact, they degrade performance for elevations between 10∘ and 35∘ .

19 of 26 American Institute of Aeronautics and Astronautics Copyright © 2013 by John W.C. Robinson, Peter Strömbäck. Published by the American Institute of Aeronautics and Astronautics, Inc., with permission.

IV.B.5.

Impact of prediction horizon on miss distance

Downloaded by John Robinson on September 5, 2013 | http://arc.aiaa.org | DOI: 10.2514/6.2013-5114

As stated in Section III.B, a path following guidance can be obtained by reinterpreting the terminal time 𝑡𝑇 as a prediction horizon. In figure 11 the effect of using shorter prediction horizons than the nominal terminal time is illustrated. For the guidance strategy developed here it can be seen that short prediction horizons are insufficient for decreasing dispersion and that the advantage of using increasingly longer prediction horizons vanishes for horizons longer than approximately 15s. It is quite possible, however, that this behavior may change if more accurate approximations are used for the higher order guidance terms than those suggested in section II.E.2 and II.E.3. Also, the poor behavior for short prediction times can probably be improved with a better path following algorithm than the one that results from our essentially predictive guidance strategy.

(a) Mean value of miss distance.

(b) Standard deviation of miss distance.

Figure 11: Miss distance as a function of prediction horizon. The projectile is launched at 30∘ elevation with initial 1𝜎-perturbation of velocity, elevation, and azimuth of 3m/s, 0.1∘ and 0.1∘ , respectively

V.

Conclusion

We have developed and analyzed a predictive guidance strategy that replaces the need for elaborate (e.g. 6DOF) on-line trajectory calculations with one high fidelity calculation done off-line combined with simple perturbation based corrections done in-flight. This is made possible by the fact that for ballistic trajectories of objects with high ballistic coefficient, such as artillery shells, the deviation in terminal state as a result of an initial perturbation is not very different from the perturbation for the corresponding vacuum trajectory (or drag-free) case. In simulations using a realistic model of a course corrected 155mm artillery shell with 2D canard control we showed that the proposed scheme was able to significantly reduce the dispersion compared to an uncorrected round, for various forms of common errors in initial values. In fact, circular error probable (CEP) values in the order of meters could be obtained. The control saturation limit was set to 80N, which can be realized with nose mounted canards of moderate size. Moreover, increasing the limit beyond 80N did not significantly improve the results. It was also shown that, for the case of shorter prediction horizon than the time-to-go, horizon lengths longer than about 15s did not improve the results. This can presumable change if better approximations to the perturbation terms are used. Indeed, including the higher order corrections in the perturbation terms (“non-vacuum trajectory corrections”) only improved the overall results slightly which indicates that the truncated expressions used in the simulations should be replaced with better approximations (using e.g. the functional equations developed). Finally, it was observed that almost all the remaining error after correction was due to muzzle velocity errors suggesting that some form of drag brake (possibly in the form of vernier correctors) could improve the 2D concept studied here.

20 of 26 American Institute of Aeronautics and Astronautics Copyright © 2013 by John W.C. Robinson, Peter Strömbäck. Published by the American Institute of Aeronautics and Astronautics, Inc., with permission.

Appendix VI. VI.A.

Approximate expressions for the acceleration evolution

Functional equations

Downloaded by John Robinson on September 5, 2013 | http://arc.aiaa.org | DOI: 10.2514/6.2013-5114

¯ in Based on the approximations introduced in section II.E.1 we shall now derive functional equations for 𝑎 ¯ that can be used in (10) and (11). (3) to be used to derive approximations for 𝑎 From the definitions (3) and (8), the drag-only approximation (14), the chain rule for differentiation and (10), (11) we obtain (recall the notation in (13)) (︀ )︀ (︀ )︀ 𝜕𝑎 Φ((𝑥, 𝑣 0 ), 𝑡), 𝑡 ⃒⃒ 𝜕𝑎𝑑 Φ((𝑥, 𝑣 0 ), 𝑡) ⃒⃒ ¯ ((𝑥, 𝑣 0 ), 𝑡) ⃒⃒ 𝜕𝑎 = = = 𝑥=𝑥0 𝑥=𝑥0 𝑥=𝑥0 𝜕𝑥 𝜕𝑥 𝜕𝑥 (︀ )︀ 𝜕𝑎𝑑 (˜ 𝑥, [Φ((𝑥0 , 𝑣 0 ), 𝑡)]4:6 ) ⃒⃒ 𝜕[Φ((𝑥, 𝑣 0 ), 𝑡)]1:3 ⃒⃒ ˜ 𝑥 =[Φ((𝑥 ,𝑣 ),𝑡)] 𝑥=𝑥0 0 0 1:3 ˜ 𝜕𝑥 𝜕𝑥 (︀ )︀ ˜) ⃒ 𝜕𝑎𝑑 ([Φ((𝑥0 , 𝑣 0 ), 𝑡)]1:3 , 𝑣 𝜕[Φ((𝑥, 𝑣 0 ), 𝑡)]4:6 ⃒⃒ ⃒ + = ˜ 𝑣 =[Φ((𝑥 ,𝑣 ),𝑡)] 𝑥=𝑥0 0 0 4:6 ˜ 𝜕𝑣 𝜕𝑥 (︀ )︀𝑇 𝜕[Φ((𝑥, 𝑣 0 ), 𝑡)]1:3 ⃒ 𝑆𝐶𝐷 𝑉 (𝑡) ⃒ − 𝑣(𝑡) ∇𝜌(𝑥(𝑡)) 𝑥=𝑥0 2𝑚 𝜕𝑥 (︀ )︀ 𝜌(𝑥(𝑡))𝑆𝐶𝐷 𝑉 (𝑡) 𝜕[Φ((𝑥, 𝑣 0 ), 𝑡)]4:6 ⃒⃒ 1 − = 𝐼 3×3 + 𝑣(𝑡)𝑣(𝑡)𝑇 2 𝑥=𝑥0 2𝑚 𝑉 (𝑡) 𝜕𝑥 ∫︁ 𝑡 )︀ (︀ )︀𝑇 (︀ ¯ ((𝑥, 𝑣 0 ), 𝑠) ⃒⃒ 𝑆𝐶𝐷 𝑉 (𝑡) 𝜕𝑎 − 𝑑𝑠 𝑣(𝑡) ∇𝜌(𝑥(𝑡)) 𝐼 3×3 + (𝑡 − 𝑠) 𝑥=𝑥0 2𝑚 𝜕𝑥 0 ∫︁ )︀ 𝑡 𝜕 𝑎 ¯ ((𝑥, 𝑣 0 ), 𝑠) ⃒⃒ 𝜌(𝑥(𝑡))𝑆𝐶𝐷 𝑉 (𝑡) (︀ 1 𝑇 𝑑𝑠, − 𝐼 3×3 + 𝑣(𝑡)𝑣(𝑡) 𝑥=𝑥0 2𝑚 𝑉 (𝑡)2 𝜕𝑥 0

𝑡 ∈ ℐ, (45)

and (︀ )︀ (︀ )︀ 𝜕𝑎 Φ((𝑥0 , 𝑣), 𝑡), 𝑡 ⃒⃒ 𝜕𝑎𝑑 Φ((𝑥0 , 𝑣), 𝑡) ⃒⃒ ¯ ((𝑥0 , 𝑣), 𝑡) ⃒⃒ 𝜕𝑎 = = = 𝑣=𝑣 0 𝑣=𝑣 0 𝑣=𝑣 0 𝜕𝑣 𝜕𝑣 𝜕𝑣 (︀ )︀ 𝜕𝑎𝑑 (˜ 𝑥, [Φ((𝑥0 , 𝑣 0 ), 𝑡)]4:6 ) ⃒⃒ 𝜕[Φ((𝑥0 , 𝑣), 𝑡)]1:3 ⃒⃒ ˜ =[Φ((𝑥0 ,𝑣 0 ),𝑡)]1:3 𝑥 𝑣=𝑣 0 ˜ 𝜕𝑥 𝜕𝑣 (︀ )︀ ⃒ ˜) ⃒ 𝜕𝑎𝑑 ([Φ((𝑥0 , 𝑣 0 ), 𝑡)]1:3 , 𝑣 𝜕[Φ((𝑥0 , 𝑣), 𝑡)]4:6 ⃒ ⃒ = ˜ =[Φ((𝑥0 ,𝑣 0 ),𝑡)]4:6 𝑣 𝑣=𝑣 0 ˜ 𝜕𝑣 𝜕𝑣 (︀ )︀𝑇 𝜕[Φ((𝑥0 , 𝑣), 𝑡)]1:3 ⃒ 𝑆𝐶𝐷 𝑉 (𝑡) ⃒ − 𝑣(𝑡) ∇𝜌(𝑥(𝑡)) 𝑣=𝑣 0 2𝑚 𝜕𝑣 )︀ 𝜕[Φ((𝑥0 , 𝑣), 𝑡)]4:6 ⃒ 𝜌(𝑥(𝑡))𝑆𝐶𝐷 𝑉 (𝑡) (︀ 1 ⃒ − 𝐼 3×3 + 𝑣(𝑡)𝑣(𝑡)𝑇 = 𝑣=𝑣 0 2𝑚 𝑉 (𝑡)2 𝜕𝑣 ∫︁ 𝑡 (︀ )︀𝑇 (︀ )︀ ¯ ((𝑥0 , 𝑣), 𝑠) ⃒⃒ 𝑆𝐶𝐷 𝑉 (𝑡) 𝜕𝑎 − 𝑣(𝑡) ∇𝜌(𝑥(𝑡)) 𝑡𝐼 3×3 + (𝑡 − 𝑠) 𝑑𝑠 𝑣=𝑣 0 2𝑚 𝜕𝑣 0 ∫︁ 𝑡 )︀(︀ )︀ ¯ 𝜌(𝑥(𝑡))𝑆𝐶𝐷 𝑉 (𝑡) (︀ 1 𝜕 𝑎 ((𝑥0 , 𝑣), 𝑠) ⃒⃒ − 𝐼 3×3 + 𝑣(𝑡)𝑣(𝑡)𝑇 𝐼 3×3 + 𝑑𝑠 , 2 𝑣=𝑣 0 2𝑚 𝑉 (𝑡) 𝜕𝑣 0

𝑡 ∈ ℐ. (46)

The relations (45), (46) represent two functional equations for the matrix functions 𝑡 ↦→

¯ ((𝑥, 𝑣 0 ), 𝑡) ⃒⃒ 𝜕𝑎 𝑥=𝑥0 𝜕𝑥

and

(47)

¯ ((𝑥0 , 𝑣), 𝑡) ⃒⃒ 𝜕𝑎 , (48) 𝑣=𝑣 0 𝜕𝑣 and different forms of approximations for the solutions to the functional equations (45), (46) shall now be derived. We note that the functions in (45) and (46) tend to zero as the ballistic coefficient tends to infinity. 𝑡 ↦→

21 of 26 American Institute of Aeronautics and Astronautics Copyright © 2013 by John W.C. Robinson, Peter Strömbäck. Published by the American Institute of Aeronautics and Astronautics, Inc., with permission.

VI.A.1.

Affine operators

If we fix a 𝑇 > 0 and define the bounded affine operator 𝐿𝑥 : 𝐶([0, 𝑇 ])3×3 → 𝐶([0, 𝑇 ])3×3 by 𝐿𝑥 𝐹 (𝑡) = −

∫︁ 𝑡 (︀ )︀𝑇 (︀ )︀ 𝑆𝐶𝐷 𝑉 (𝑡) 𝑣(𝑡) ∇𝜌(𝑥(𝑡)) 𝐼 3×3 + (𝑡 − 𝑠)𝐹 (𝑠) 𝑑𝑠 2𝑚 0 ∫︁ )︀ 𝑡 𝜌(𝑥(𝑡))𝑆𝐶𝐷 𝑉 (𝑡) (︀ 1 𝑇 𝑣(𝑡)𝑣(𝑡) 𝐹 (𝑠) 𝑑𝑠, − 𝐼 3×3 + 2𝑚 𝑉 (𝑡)2 0

𝑡 ∈ [0, 𝑇 ], (49)

we see from (45) that the matrix function in (47) is a solution to the equation 𝐹 = 𝐿𝑥 𝐹

(50)

Downloaded by John Robinson on September 5, 2013 | http://arc.aiaa.org | DOI: 10.2514/6.2013-5114

i.e. the function in (47) is a fix point to the operator 𝐿𝑥 . Analogous arguments can be applied to the operator 𝐿𝑣 : 𝐶([0, 𝑇 ]3×3 ) → 𝐶([0, 𝑇 ]3×3 ) defined by 𝐿𝑣 𝐹 (𝑡) = −

∫︁ 𝑡 (︀ )︀𝑇 (︀ )︀ 𝑆𝐶𝐷 𝑉 (𝑡) 𝑣(𝑡) ∇𝜌(𝑥(𝑡)) 𝑡𝐼 3×3 + (𝑡 − 𝑠)𝐹 (𝑠) 𝑑𝑠 2𝑚 0 ∫︁ 𝑡 )︀ )︀(︀ 𝜌(𝑥(𝑡))𝑆𝐶𝐷 𝑉 (𝑡) (︀ 1 𝑇 𝐹 (𝑠) 𝑑𝑠 , 𝐼 3×3 + 𝑣(𝑡)𝑣(𝑡) 𝐼 + − 3×3 2 2𝑚 𝑉 (𝑡) 0

𝑡 ∈ [0, 𝑇 ]. (51)

Here we see that the function in (48) is a solution to the equation 𝐹 = 𝐿𝑣 𝐹 ,

(52)

i.e. the function in (48) is a fix point to the operator 𝐿𝑣 . VI.A.2.

Short-time properties

From the right hand sides of (45) and (46) we straightforwardly get two estimates, valid for small 𝑡, of the functions in (47) and (48) as

and

)︀𝑇 ¯ ((𝑥, 𝑣 0 ), 𝑡) ⃒⃒ 𝜕𝑎 𝑆𝐶𝐷 𝑉0 (︀ =− 𝑣 0 ∇𝜌(𝑥0 ) + 𝒪(𝑡), 𝑥=𝑥 0 𝜕𝑥 2𝑚

(53)

)︀ ¯ ((𝑥0 , 𝑣), 𝑡) ⃒⃒ 𝜌(𝑥0 )𝑆𝐶𝐷 𝑉0 (︀ 1 𝜕𝑎 =− 𝐼 3×3 + 2 𝑣 0 𝑣 𝑇0 + 𝒪(𝑡), 𝑣=𝑣 0 𝜕𝑣 2𝑚 𝑉0

(54)

where 𝒪 is big O (ordo) and 𝑉0 = ‖𝑣 0 ‖. The two simple approximations in (53) and (54) shall be used below as a start for deriving approximations of the functions in (47) and (48). VI.B.

Iterative approximation

The operators 𝐿𝑥 , 𝐿𝑣 in (49) and (51), respectively, will be contractive under several different conditions, the most obvious being when 𝑇 is sufficiently smalls since ∫︁ (︀ )︀ (︀ )︀𝑇 𝑡 𝑆𝐶𝐷 𝑉 (𝑡) 𝐿 𝐹 1 (𝑡) − 𝐹 2 (𝑡) = − 𝑣(𝑡) ∇𝜌(𝑥(𝑡)) (𝑡 − 𝑠)(𝐹 1 (𝑠) − 𝐹 2 (𝑠)) 𝑑𝑠 2𝑚 0 ∫︁ )︀ 𝑡 𝜌(𝑥)𝑆𝐶𝐷 𝑉 (𝑡) (︀ 1 𝑇 − 𝐼 3×3 + 𝑣(𝑡)𝑣(𝑡) (𝐹 1 (𝑠) − 𝐹 2 (𝑠)) 𝑑𝑠 (55) 2𝑚 𝑉 (𝑡)2 0 for both 𝐿 = 𝐿𝑥 and 𝐿 = 𝐿𝑣 , and from this an estimate for the coefficient of contraction can easily be obtained (see below). s The

operators 𝐿𝑥 , 𝐿𝑣 are also contractive when the ballistic coefficient is sufficiently high.

22 of 26 American Institute of Aeronautics and Astronautics Copyright © 2013 by John W.C. Robinson, Peter Strömbäck. Published by the American Institute of Aeronautics and Astronautics, Inc., with permission.

When the operators 𝐿𝑥 , 𝐿𝑣 are contractive, the two equations (50) and (52) have unique solutions and these can be found by iterative approximationt using the scheme 𝐹 𝑛+1 𝐹0

= {︃ 𝐿𝐹 , 𝑛 = 0, 1, 2, . . . ⃒ ¯ ((𝑥,𝑣 0 ),0) ⃒ 𝜕𝑎 if 𝜕𝑥 ⃒𝑥=𝑥0 = ¯ ((𝑥0 ,𝑣),0) ⃒ 𝜕𝑎 if 𝜕𝑣 𝑣=𝑣 0

𝐿 = 𝐿𝑥 ,

(56)

𝐿 = 𝐿𝑣 ,

provided we can guarantee that the iterates fall within the domain of attraction of the fix point (i.e. the domain where 𝐿 is contractive). From the size estimates in (18)–(20) and the relation (55) we straightforwardly get the estimate (︀ )︀ ‖𝐿 𝐹 1 (𝑡) − 𝐹 2 (𝑡) ‖ ≤ 0.03(𝑡 + 10−2 𝑡2 )‖𝐹 1 (𝑡) − 𝐹 2 (𝑡)‖ and therefore 𝐿𝑥 , 𝐿𝑣 are contractive at least for 𝑡 ≤ 25s for the standard 155mm shell firing case. Moreover, for 𝑡 ≤ 10s the Lipschitz constant is less than 1/3 so the convergence in the iterations in (56) is fast.

Downloaded by John Robinson on September 5, 2013 | http://arc.aiaa.org | DOI: 10.2514/6.2013-5114

VI.C.

First order approximation

If we use the short-time approximations obtained from (53), (54) (by neglecting higher order terms) as zeroth order approximations in the iteration scheme (56) we get as first order approximation to the fix point for 𝐿𝑥 the function (recall the notation (13)) 𝐹 𝑥,1 (𝑡) = −

(︀ )︀𝑇 (︀ )︀𝑇 )︀ 𝑡2 𝑆𝐶𝐷 𝑉0 (︀ 𝑆𝐶𝐷 𝑉 (𝑡) 𝑣(𝑡) ∇𝜌(𝑥(𝑡)) 𝐼 3×3 − 𝑣 0 ∇𝜌(𝑥0 ) 2𝑚 2 2𝑚 2 )︀ (︀ )︀𝑇 𝜌(𝑥(𝑡))𝑆 2 𝐶𝐷 𝑉0 𝑉 (𝑡) (︀ 1 +𝑡 𝐼 3×3 + 𝑣(𝑡)𝑣(𝑡)𝑇 𝑣 0 ∇𝜌(𝑥0 ) , 2 2 4𝑚 𝑉 (𝑡)

𝑡 ∈ ℐ, (57)

and as first order approximation to the fix point for 𝐿𝑣 we get (︀ )︀𝑇 (︁ )︀)︁ 𝑆𝐶𝐷 𝑉 (𝑡) 1 𝑡2 𝜌(𝑥0 )𝑆𝐶𝐷 𝑉0 (︀ 𝑣(𝑡) ∇𝜌(𝑥(𝑡)) 𝐼 3×3 + 2 𝑣 0 𝑣 𝑇0 𝑡𝐼 − 2𝑚 2 2𝑚 𝑉0 (︁ )︀ )︀)︁ 𝜌(𝑥(𝑡))𝑆𝐶𝐷 𝑉 (𝑡) (︀ 1 𝜌(𝑥0 )𝑆𝐶𝐷 𝑉0 (︀ 1 𝑇 𝑇 − 𝐼 3×3 + 𝑣(𝑡)𝑣(𝑡) 𝐼 − 𝑡 𝐼 + 𝑣 𝑣 , 3×3 3×3 0 0 2 2𝑚 𝑉 (𝑡)2 2𝑚 𝑉0

𝐹 𝑣,1 (𝑡) = −

VI.C.1.

𝑡 ∈ ℐ. (58)

Using left endpoint values

For off line computation the operators 𝐿𝑥 , 𝐿𝑣 in (49), (51) and associated approximations in (57), (58) should present few problemsu but for real time computations by an on board guidance computer the complexity of the involved operations may be too large. One way to simplify the computations is to redefine the operators 𝐿𝑥 , 𝐿𝑣 so that they only contain left endpoint values (𝑥0 , 𝑣 0 ), and again start the iterations in (56) with the constant short-time approximations obtained from (53), (54) as zeroth order approximation. With this approach, the right hand sides of (57), (58) and the approximate solutions to the functional equations (50) and (52) will only contain left endpoint values (𝑥0 , 𝑣 0 ) and be polynomial in time. The integrals in (47) and (48) will then be trivial to compute. If this approach is applied and only terms of order 𝑡2 or lower are retained, and the result is inserted in the integrals on the right hand side in (10) and (11), the approximations in (28)–(31) are obtained. VI.D.

Alternative methods

Another method to obtain approximations to the integrals on the right hand side in (10) and (11) is to employ some or all terms of the approximations in (57) and (58) together with a numerical quadrature scheme. For instance, a simple three point Gauss-Lobatto scheme may be used, where the values at the rightmost node point is taken from the nominal trajectory. Several other approaches are possible. t One could equally well combine (10), (11) and (45), (46) to set up a contraction mapping iteration scheme directly for ¯. 𝜕Φ((𝑥0 , 𝑣 0 ), 𝑡)/𝜕𝑥0 and 𝜕Φ((𝑥0 , 𝑣 0 ), 𝑡)/𝜕𝑥0 , but it is here more convenient to work with 𝑎 u An alternative method of computing the functions in (47) and (48) off line is to compute them directly from e.g. a 6DOF calculation.

23 of 26 American Institute of Aeronautics and Astronautics Copyright © 2013 by John W.C. Robinson, Peter Strömbäck. Published by the American Institute of Aeronautics and Astronautics, Inc., with permission.

VII.

Guidance command filter

Downloaded by John Robinson on September 5, 2013 | http://arc.aiaa.org | DOI: 10.2514/6.2013-5114

The guidance commands are transformed from the fire control frame 𝐹 to the non spinning frame 𝑁𝑠 (in which the autopilot works) and then filtered through a second order linear filter with a rate limiter before being sent to the (static) autopilot (see below). The object of the linear filter is to smooth out the guidance commands when they are small and the rate limiter is in place to guarantee that larger changes in commands can never be passed on directly to actuators (since the precession and nutation modes are weakly damped). More precisely, each “channel” of the guidance acceleration command, i.e. the command 𝑎𝑧 or 𝑎𝑦 in (60) below for acceleration in either 𝑦 or 𝑧 direction in 𝑁𝑠 , is filtered through ]︃ ]︃ [︃ ]︃ [︃ [︃ ]︃ [︃ ˙ −𝑟ℓ ), 𝑟ℓ ) 0 0 1 min(max(𝜉, 𝜉˙ min(max(𝑎𝑐 , −𝑐ℓ ), 𝑐ℓ ), + = 𝜔𝑔2 −𝜔𝑔2 −2𝜔𝑔 𝜁𝑔 𝜉 𝜉¨ where 𝜔𝑔 , 𝜁𝑔 > 0 are the damping and natural (angular) frequency parameters, respectively, and 𝑐ℓ , 𝑟ℓ > 0 is a value and rate limit, respectively. (For open loop response characteristics of a similar shell-fuze combination, see Ref. 35.) The filtered acceleration command 𝑎𝑐,𝑓 passed to the autopilot is given by 𝑎𝑐,𝑓 = 𝜉 (i.e. the two components in 𝑎1 in (60) below are obtained from the output of the guidance filter when the guidance and control system is engaged).

VIII.

Autopilot

The autopilot described here performs only the most basic function of translating acceleration commands to steady-state values for the control forces produced by the actuators. VIII.A.

Controlled dynamics in 𝑁𝑠

The dynamics of a spinning course corrected projectile are naturally formulatedv in the non spinning frame 𝑁𝑠 . Let 𝑎 = [𝑎𝑥 , 𝑎𝑦 , 𝑎𝑧 ]𝑇 be given by 𝑎 = (𝑅(𝑁𝑠 ,𝐹 ) )𝑇 (𝑎 + 𝑔 + 𝑢),

(59)

where 𝑎, 𝑔 and 𝑢 are defined in (2), and 𝑅(𝑁𝑠 ,𝐹 ) is the rotation matrix for transformation from the non spinning frame 𝑁𝑠 to the fire control frame 𝐹 , and define 𝑎1 by 𝑎1 = [𝑎𝑧 , 𝑎𝑦 ]𝑇 .

(60)

Thus, the vector 𝑎1 holds the 𝑧 and 𝑦-components in 𝑁𝑠 of the sum of aerodynamic, gravity and possibly control induced acceleration. When the control acceleration 𝑢 is given by the predictive guidance scheme ˆ in (44). Next, let 𝑢1 ∈ R2 be given by developed in section III then 𝑢 is the command filtered version of 𝑢 𝑢1 = [𝐹𝑐,𝑧 , 𝐹𝑐,𝑦 ]𝑇 , where 𝐹𝑐,𝑧 , 𝐹𝑐,𝑦 are the 𝑧 and 𝑦-components, respectively, in 𝑁𝑠 of the control forcesw (applied near the tip of the nose). By linearizing the equations of motion in 𝑁𝑠 for the projectile (with fuze) around the value 𝑥 = 0 of the state vector 𝑥 = [𝛼, 𝛽, 𝑞, 𝑟]𝑇 , where 𝛼, 𝛽 are the angle of attack and sideslip angle in 𝑁𝑠 , and 𝑞, 𝑟 are the pitch and yaw angular rates in 𝑁𝑠 , respectively, the equations of motion can be written 𝑥˙ 𝑎1

= 𝐴(𝑞𝑑 , 𝑉, 𝑝)𝑥 + 𝑔 𝑓 (Θ, 𝑉 ) + 𝐵(𝑉 )𝑢1 , 1 = 𝐶(𝑞𝑑 , 𝑉, 𝑝)𝑥 + 𝑔 1 (Θ) + 𝑢1 , 𝑚

v The

(61) (62)

linearized dynamics21, 36 have the same form in 𝑁𝑠 as in the nonrolling frame 𝑁𝑟 since the two frames differ only by a rotation about the main axis of the projectile.31 w Note that, despite the notation, the control force vector 𝑢 in 𝑁 is only implicitly related to the control acceleration vector 𝑠 1 𝑢 in 𝐹 in the model (2). The relationship is that the sum of the three accelerations 𝑔 + 𝑎 + 𝑢 should transform as (59).

24 of 26 American Institute of Aeronautics and Astronautics Copyright © 2013 by John W.C. Robinson, Peter Strömbäck. Published by the American Institute of Aeronautics and Astronautics, Inc., with permission.

where ⎡

𝑑 𝑆𝑝𝑑 − 𝑞𝑑 𝑆 (𝐶𝐿𝛼 + 𝐶𝐷 ) − 𝑞𝑚𝑉 2 𝐶𝑁𝑝𝛼 ⎢ 𝑚𝑉𝑞𝑑 𝑆𝑝𝑑 𝑞𝑑 𝑆 𝐶 − (𝐶 ⎢ 𝑁𝑝𝛼 𝐿𝛼 + 𝐶𝐷 ) 𝑚𝑉 2 𝑚𝑉 ⎢ 𝑞𝑑 𝑆𝑑 𝑞𝑑 𝑆𝑝𝑑2 𝐴(𝑞𝑑 , 𝑉, 𝑝) = ⎢ 𝐶𝑀𝑝𝛼 (𝐵) 𝐶𝑀𝛼 (𝐵) ⎢ 𝐽𝑦𝑦 𝐽𝑦𝑦 𝑉 ⎣ 2 𝑞𝑑 𝑆𝑝𝑑 𝑆𝑑 𝐶𝑀𝑝𝛼 − 𝑞𝑑(𝐵) 𝐶𝑀𝛼 (𝐵)

𝐽𝑦𝑦 𝑉

⎤ + 𝐶𝑁𝛼˙ ) ⎥ −1 ⎥ ⎥ (𝐵) 𝐽𝑥𝑥 ⎥ , (63) −𝑝 (𝐵) ⎥ 𝐽𝑦𝑦 ⎦ 𝑞𝑑 𝑆𝑑2 (𝐶𝑀𝑞 + 𝐶𝑀𝛼˙ ) (𝐵) 𝑞𝑑 𝑆𝑑 𝑚𝑉 2 (𝐶𝑁𝑞

1 𝑞𝑑 𝑆𝑑 − 𝑚𝑉 ) 2 (𝐶𝑁𝑞 + 𝐶𝑁𝛼 ˙ 𝑞𝑑 𝑆𝑑2 (𝐵) 𝐽𝑦𝑦 𝑉

𝐽𝑦𝑦

(𝐶𝑀𝑞 + 𝐶𝑀𝛼˙ ) 𝑝

(𝐵) 𝐽𝑥𝑥 (𝐵)

𝐽𝑦𝑦

𝐽𝑦𝑦 𝑉

[︃ 𝐶(𝑞𝑑 , 𝑉, 𝑝) = 𝑉 [𝐼 2×2 , 02×2 ]𝐴(𝑞𝑑 , 𝑉, 𝑝) + 𝑉 [02×2 , 𝑆],

𝑆=

−1 0

0 1

]︃ ,

(64)

and [︃

Downloaded by John Robinson on September 5, 2013 | http://arc.aiaa.org | DOI: 10.2514/6.2013-5114

𝐵(𝑉 ) =

1 𝑚𝑉 𝐼 2×2 ℓ (𝐵) 𝑆 𝐽𝑦𝑦

]︃ ,

𝑔 𝑓 (Θ, 𝑉 ) = [

𝑔 cos(Θ), 0, 0, 0]𝑇 , 𝑉

𝑔 1 (Θ) = [𝑔 cos(Θ), 0]𝑇 .

Here, 𝑞𝑑 is the dynamic pressure, 𝑉 is the airspeed, 𝑝 is the spin rate in the body fixed frame 𝐵, 𝑚 is the (𝐵) mass, 𝐽𝑦𝑦 is the moment of inertia about the 𝑦-axis in 𝐵 and ℓ > 0 is the (forward) distance from the control force application point to CoM. The quantities 𝐶( · ) are aerodynamic coefficients (in the BRL notation, cf. Ref. 26, tab. 2.3), Θ is the Euler pitch angle (between the 𝑥-axis in 𝐵 and the horizontal plane in 𝐹 ) and 𝑔 is the constant of acceleration due to Earth’s gravity. Henceforth we shall for simplicity drop the arguments of 𝐴, 𝐶, 𝐶, 𝑔 𝑓 and 𝑔 𝑓 , 𝑔 1 . VIII.A.1.

Trimmed acceleration

From (61) and (62) we see that the problem of maintaining a given (total) acceleration 𝑎1 as in (60) can be written ]︃ [︃ ]︃ [︃ ]︃ [︃ −𝑔 𝑓 𝑥 𝐴 𝐵 . (65) = 1 𝐼 2×2 𝑢1 𝐶 𝑚 𝑎1 − 𝑔 1 Clearly, as long as the matrix on the right is nonsingular, then for every quadruple (𝑎1 , Θ, 𝑉, 𝑝) there exists a unique pair (𝑥, 𝑢1 ) such that (65) holds. Since the lower right 2 × 2 corner of the matrix on the right in (65) is nonsingular it follows, by a Schur complement argument, that the entire matrix is nonsingular if and only if 𝐴 − 𝑚𝐵𝐶 is nonsingular.x Based on this it is straightforward to show that the matrix on the right in (65) is in general nonsingular for conventional projectiles and the trim equations can then be solved. Realizing trimmed acceleration. The guidance commands are translated by the autopilot into control force values 𝐹𝑐,𝑧 , 𝐹𝑐,𝑦 that give the desired acceleration 𝑎1 at trim, i.e. at the equilibrium condition represented by the first row of (65). For 𝑢 ̸= 0 in (2) and (59) this corresponds to flight at a different state 𝑥 than at yaw of repose for the projectile.

Acknowledgments This work was supported by the Swedish Armed Forces research and technology program (Weapons and Protection).

References 1 Dornblut,

M., “Next-Generation Artillery Ammunition,” Military Technology (MILTECH), , No. 12, 2010, pp. 61–65. J. and Siewart, J., “2-D Projectile Trajectory Corrector,” January 2003, U.S. Patent No.: 6,502,786 B2. 3 Ziliani, A., Grignon, C., Trouillot, C., and Jeannin, “Diagnostic of the Behaviour of a Course-Correction Ammunition During its Correction Phase,” Proc. 19:th Intl. Symposium of Ballistics, Interlaken, Switzerland, May 2001, pp. 479–487. 4 Pettersson, T., Buretta, R., and Cook, D., “Aerodynamics and Flight Stability for a Course Corrected Artillery Round,” Proc. 23:rd Intl. Symposium on Ballistics, Tarragona, Spain, April 2007, pp. 747–754. 5 Grignon, C., Cayzac, R., and Heddadj, S., “Improvement of Artillery Projectile Accuracy,” Proc. 23:rd Intl. Symposium on Ballistics, Tarragona, Spain, April 2007, pp. 747–754. 2 Rupert,

x It is a standard fact from the theory of Schur complements that the matrix [𝐴, 𝐵; 𝐶, 𝐷] (Matlab notation) is nonsingular if and only if both 𝐷 and 𝐴 − 𝐵𝐷−1 𝐶 are nonsingular (this follows e.g. from the Schur determinant formula, see e.g.37 ).

25 of 26 American Institute of Aeronautics and Astronautics Copyright © 2013 by John W.C. Robinson, Peter Strömbäck. Published by the American Institute of Aeronautics and Astronautics, Inc., with permission.

Downloaded by John Robinson on September 5, 2013 | http://arc.aiaa.org | DOI: 10.2514/6.2013-5114

6 Jermey, C., “Wind Tunnel Tests of a Spinning 105mm Artillery Shell Model With Canard Control Surfaces,” Tech. Rep. WSRL-0122-TR, Weapons Systems Research Laboratory, Dept. Defence, Defence Science and Technology Organisation (DSTO), Adelaide, South Australia, September 1979. 7 Regan, F. and Smith, J., “Aeroballistics of a Terminally Corrected Spinning Projectile (TCSP),” Journal of Spacecraft, Vol. 12, No. 12, December 1975, pp. 733–738. 8 Clancy, J., Bybee, T., and Friedrich, W., “Fixed Canard 2-D Guidance of Artillery Projectiles,” January 2006, U.S. Patent No.: 6,981,672 B2. 9 Sebestyen, G., Sinclair, R., Smith, J., Sands, T., and Nussdorfer, T., “Canard Control Assembly for a Projectile,” April 1985, U.S. Patent No.: 4,512,537. 10 Ollerenshaw, D. and Costello, M., “Simplified Projectile Swerve Solution for General Control Inputs,” Journal of Guidance, Control and Dynamics, Vol. 31, No. 5, Sept.-Oct. 2008, pp. 1259–1265. 11 Fresconi, F. and Plostins, P., “Control Mechanism Strategies for Spin-Stabilized Projectiles,” Tech. Rep. ARL-TR-4611, US Army Research Laboratory, September 2008. 12 Bakken, R., “Modellering av roterende legeme (Modeling of rotating body),” Tech. Rep. FFI-rapport 2009/01503, Norwegian Defence Research Establishment (FFI), September 2009, (In Norwegian). 13 Lloyd, K. and Brown, D., “Instability of Spinning Projectiles During Terminal Guidance,” Journal of Guidance and Control, Vol. 2, No. 1, Jan-Feb 1979, pp. 65–70. 14 Murphy, C., “Instability of Controlled Projectiles in Ascending or Descending Flight,” Journal of Guidance and Control, Vol. 4, No. 1, Jan-Feb 1981, pp. 66–69, Also as AIAA paper 79-1669. 15 Fresconi, F., Brown, T., Celmins, I., DeSpirito, J., Ilg, M., and Maley, J., “Very Affordable Precision Projectile System and Flight Experiments,” Proc. 27:th Army Science Conference, Orlando, FL, 2010, Paper DO-04. 16 Andersson, K., Bartelson, N., and Bondesson, S., “Common Geometry Smart Projectiles, Early Studies and Tests of A Concept,” Proc. Eighth Intl. Symposium on Ballistics, Orlando, FL, October 1984, pp. IV 57–68. 17 Gagnon, E. and Lauzon, M., “Course Correction Fuze Concept Analysis for In-Service 155 mm Spin-Stabilized Gunnery Projectiles,” Proc. AIAA Guidance, Navigation and Control Conference and Exhibit, Honolulu, HI, August 2008, AIAA paper 2008-6997. 18 Ollerenshaw, D. and Costello, M., “Model Predictive Control of a Direct Fire Projectile Equipped with Canards,” Proc. AIAA Guidance, Navigation, and Control Conference and Exhibit, San Francisco, CA, August 2005, AIAA paper 2005-5818. 19 Siouris, G., Missile Guidance and Control Systems, Springer-Verlag, New York, NY, 2004. 20 Gagnon, E. and Lauzon, M., “Maneuverability Analysis of the Conventional 155 mm Gunnery Projectile,” Proc. AIAA Guidance, Navigation and Control Conference and Exhibit, Hilton Head, SC, August 2007, AIAA paper 2007-6784. 21 Calise, A. and El-Shirbiny, H., “An Analysis of Aerodynamic Control for Direct Fire Spinning Projectiles,” Proc. AIAA Guidance, Navigation and Control Conference and Exhibit, Montreal, Canada, August 2001, AIAA paper 2001-4217. 22 Theodoulis, S. and Wernert, P., “Flight Control for a Class of 155 mm Spin-Stabilized with Course Correction Fuze (CCF),” Proc. AIAA Guidance, Navigation and Control Conference and Exhibit, Portland, OR, August 2011, AIAA paper 2011-6247. 23 Theodoulis, S., Grassmann, V., Wernert, P., Dritsas, L., Kitsios, I., and Tzes, A., “Guidance and Control Design for a Class of Spin-Stabilized Fin-Controlled Projectiles,” J. Guidance, Control and Navigation, Vol. 36, No. 2, March-April 2013, pp. 517–531. 24 Salman, M. and Chang, B., “Active Coning Compensation for Control of Spinning Flying Vehicles,” Proc. IEEE Intl. Conference on Control Applications, Part of IEEE Multi-Conference on Systems and Control, Yokohama, Japan, September 2010. 25 Fresconi, F., Cooper, G., and Costello, M., “Practical Assessment of Real-Time Impact Point Estimators for Smart Weapons,” J. Aerospace Engineering, Vol. 24, No. 1, pp. 1–11. 26 McCoy, R., Modern Exterior Ballistics: The Launch and Flight Dynamics of Symmetric Projectiles, Schiffer Publishing, Atglen, PA, 1999, See also the errata by D.G. Miller et al. available on internet. 27 Alekseev, V., “An Estimate for the Perturbations of the Solutions of Ordinary Differential Equations,” Vestn. Moskov. Univ. Ser. I Mat. Meh., , No. 2, 1961, pp. 28–36, (In Russian). 28 Brauer, F., “Perturbations of Nonlinear Systems of Differential Equations. II,” J. Math. Anal. Appl., Vol. 17, 1967, pp. 418–434. 29 Qu, Z., Robust Control of Nonlinear Uncertain Systems, Nonlinear Science, Wiley, New York, NY, 1998. 30 Johnson, D., Minimum Distance Quesries for Haptic Rendering, Ph.D. thesis, Univ. Utah, May 2005. 31 Stilley, G. and Alford, R., “Axis System Considerations for Guidance and Stability of Spinning Projectiles,” Proc. AIAA 22nd Aerospace Sciences Meeting, Reno, NV, January 1984, AIAA-paper 84-0327. 32 Khalil, M., Abdalla, H., and Kamal, O., “Dispersion Analysis for Spinning Artillery Projectile,” Proc. 13:th Intl. Conference on Aerospace Science and Aviation Technology, May 2009, ASAT-13-FM-03. 33 Altshuler, B., “Report on the Firing Test to Determine the Limits of Accuracy Life of 4.7” A.A. Gun T2,” Tech. Rep. BRL R 254, Ballistics Research Laboratory, Aberdeen Proving Ground, MD, September 1941. 34 Roberts, N., “Ballistic Analysis of Firing Table Data For 155mm, M825 Smoke Projectile,” Tech. Rep. BRL-MR-3865, Ballistics Research Laboratory, Aberdeen Proving Ground, MD, September 1990. 35 Wernert, P., Leopold, F., Bidino, D., Junker, J., Lehmann, L., B¨ ar, K., and Reindler, A., “Wind tunnel tests and openloop trajectory simulations for a 155 mm canards guided spin stabilized projectile,” Proc. AIAA Atmoshperic Flight Mechanics Conference and Exhibit, Honolulu, HI, August 2008, AIAA paper 2008-6881. 36 Mracek, C., Stafford, M., and Unger, M., “Control of Spinning Symmetric Airframes,” Proc. AIAA Missile Sciences Conference, Monterey, CA, November 2006. 37 Cottle, R., “Manifestations of the Schur Complement,” Linear Algebra and Its Applications, Vol. 8, 1974, pp. 189–211.

26 of 26 American Institute of Aeronautics and Astronautics Copyright © 2013 by John W.C. Robinson, Peter Strömbäck. Published by the American Institute of Aeronautics and Astronautics, Inc., with permission.