AIAA 2015-0087 AIAA SciTech 5-9 January 2015, Kissimmee, Florida AIAA Guidance, Navigation, and Control Conference

Velocity To Be Gained Guidance for a Generic 2D Course Correcting Fuze John W.C. Robinson* and Peter Str¨omb¨ack†

Downloaded by John Robinson on January 26, 2015 | http://arc.aiaa.org | DOI: 10.2514/6.2015-0087

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

A guidance strategy for generic 2D course correcting fuzes is developed using a perturbation based trajectory prediction method together with a velocity to be gained (VTG) control formulation. The open loop dynamics for the VTG that results in the atmospheric case is a slight modification of the vacuum dynamics and similar approaches for control can be applied. We formulate and solve a time varying linear quadratic (LQ) optimal control problem to obtain a guidance strategy that makes efficient use of controls over the flight. The resulting guidance strategy requires only low complexity updates in-flight. Using simulations of a standard 155mm artillery shell with canard control in varying firing scenarios it is shown that excellent dispersion reduction can be obtained with only modest requirements on control authority.

I.

Introduction

Long range cannons for indirect fire is a central component in most field artillery systems. With a typical 155mm artillery system the achievable range is well over 20km, even without using base bleed or rocket assist for range enhancement. At these ranges the dispersion in impact point of standard unguided rounds can be several times the lethal radius (about 50m) of a typical high explosive fragmentation type shell, at least in the downrange direction. Since this dispersion figure has remained essentially unchanged for a long time it appears that the only realistic route to significantly reducing dispersion is by introducing guidance and control capabilities in the projectiles. In the last few decades several variants of cannon fired projectiles with such capabilities have been developed and fielded. The first examples of this development were specially designed projectiles, offering true precision weapon capabilities (high probability of target destruction with a single round). However, the high cost of these rounds has somewhat hindered their wide deployment and more cost effective solutions have been sought. One way to realize more cost effective solutions is via “guidance kits”. A guidance kit is a plug-in replacement for the ordinary fuze which provides course correction capabilities by employing e.g. canards or drag brakes. The fact that cannon based artillery systems typically have high nominal accuracy and precision, and that true precision weapon capability might not always be needed, means that a solution with a fuze that is capable of (small) course corrections can often represent the best trade off between performance and cost. Moreover, course correcting fuzes can be fitted on varying portions of the large existing stockpile of standard shells, thereby offering a gradual “on demand” path to increased precision. I.A.

Controlling dispersion with course correcting fuzes

For unguided artillery projectiles, the major contributing factors to dispersion 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. In applications where high accuracy and precision is needed in both downrange and crossrange directions, fuzes * Deputy † Senior

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

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

with true 2D correction capability must be used. The standard way of achieving good 2D course correction is via nose mounted canards which can generate normal (longitudinal and lateral) control forces, with a controllable magnitude. The idea of using nose mounted canards for generation of 2D control forces on course corrected projectiles (e.g. via a fuze kit) is not new. A design based on a combination of movable and fixed canards was proposeda and analyzed in the 1970’s by Regan and Smith.2 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 spin while the remaining pair of canards can be deflected to produce a net normal force in a direction given by the fuze rotation angle. Variants of this idea have subsequently been implemented in existing course correcting fuzes. We shall here consider a generic variant of the 2D canard control concept, where the basic assumption is that the normal force exerted by the canards at the nose can be directed at any direction at any time and continuously adjusted,3 thus neglecting to model in detail how this is actually accomplished.

Downloaded by John Robinson on January 26, 2015 | http://arc.aiaa.org | DOI: 10.2514/6.2015-0087

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 this type of control for 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.4 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.5, 4 Furthermore, a too large normal force near the tip of the projectile can easily lead to lossb of stability.7, 8 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 if stability is to be guaranteed by the natural stability of the projectile. Restricting the control to course corrections translates to small control effort and for this predictive guidance is particularly well suited. 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.9 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. Guidance based on (somewhat sophisticated) prediction thus offers the possibility of low control effort but is susceptible to modeling errors and changing environmental conditions. Therefore, in reality at least a few such computations with accompanying corrections of the path are needed during the flight with predictive guidance. As a consequence, the inherently open loop predictive guidance strategy will in practice lead to a closed loop one. In the extreme case a closed loop guidance strategy can be formed where the predictions are recomputed repeatedly. This includes the velocity to be gained guidance we develop below where the the guidance loop is described by a dynamical relation (a differential equation) and the error quantity that drives the guidance loop is based on a prediction of the terminal state which is recomputed repeatedly (in a continuous time fashion). I.B.1.

Related work

Both predictive guidance and path following guidance for spinning projectiles have received some interest in the literature. Gagnon and Lauzon10, 11 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 found11 that drag brakes could effectively correct muzzle velocity errors but canards were superior for aiming errors and crosswind errors. Calise and El-Shirbiny12 studied predictive guidance for direct fire spinning projectiles and noted 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. Gross and Costello13 developed a predictive guidance and control algorithm with protection against control induced instability. The algorithm was based on a model predictive control formulation using a nonlinear six degree of freedom model for the projectile a An

early aerodynamic investigation of a 105mm projectile with movable canards can be found in Ref. 1. 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.6 b Another

2 of 24 American Institute of Aeronautics and Astronautics

Downloaded by John Robinson on January 26, 2015 | http://arc.aiaa.org | DOI: 10.2514/6.2015-0087

dynamics. Ollerenshaw and Costello9 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 inflight using so called projectile linear theory. It was found that prediction horizon had a large effect on the results, with a longer horizon leading to reduced dispersion. The stability and disturbance rejection for 2D canard controlled spinning projectiles can be improved by incorporation of a dynamics synthesizing autopilot and several such designs have been presented recently.14, 15, 16 It has been shown that good tracking performance can be achieved e.g. 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. 17. The velocity to be gained concept we shall use below was originally developed for ballistic missile and spaceflight applications in the 1950’s by Laning and Battin at MIT’s instrumentation laboratory.18, 19, 20, 21 Velocity to be gained has remained a useful tool to develop guidance solutions and recent applications include high satellite orbit insertion22, 23 and ballistic missile defense.24 The idea of applying perturbation expansions to develop guidance strategies also goes back to the era when the velocity to be gained concept was conceived.20 I.B.2.

Contribution

In the current work we extend some of the perturbation guidance results for the atmospheric case that were developed specifically for spinning projectile guidance applications in Ref. 25. Using these results as a basis, we propose a version of velocity to be gained guidance which is adapted to course correcting fuzes with generic 2D control. The derivation includes a fairly detailed treatment of the assumptions and approximations that can be employed in the atmospheric (constant gravity) case. For the closed loop guidance we formulate and solve a linear quadratic control problem which is designed to penalize errors in the guidance loop in such a way that control action is performed early in the flight, when it provides large “leverage”. The resulting algorithm requires only low complexity on-line computations during flight and is highly effective in reducing the dispersion. Moreover, it is capable of doing so with very modest actuator requirements due to the fact that the control effort is distributed evenly over the flight phase where the guidance is active. I.C.

Notation

Notation is standard and we shall largely omit units. In general, we consider points in R𝑛 as column vectors. The norm ‖ · ‖ is always the 2-norm (for a vector, and the induced 2-norm for a matrix). Points in R𝑛 × R𝑛 will sometimes be denoted as (𝑥, 𝑦) and sometimes as [𝑥𝑇 , 𝑦 𝑇 ]𝑇 where superscript 𝑇 denotes transpose. If 𝑀 is a matrix or a vector, the notation [𝑀 ]𝑗:𝑘 denotes the 𝑗th through 𝑘th rows (inclusive). The plane perpendicular to 𝑥 ∈ R3 is marked [𝑥]⊥ . For a vector 𝜉 depending on another vector 𝜂 the Jacobian matrix with entries 𝜕𝜉 𝑖 /𝜕𝜂 𝑖 is denoted 𝜕𝜉/𝜕𝜂. The identity matrix in R3×3 is marked 𝐼 3×3 and 0 is a zero vector. The symbol 𝑜 is (Landau’s) small o (ordo) and 𝒪 is big O. A function (scalar, vector or matrix) which is 𝑘 ≥ 0 times continuously differentiable (in any element and with respect to any argument) is said to be 𝐶 𝑘 . I.D.

Outline

In the next section we introduce some notation and describe our trajectory model and its associated approximations. After this, the proposed guidance strategy is developed in section III, beginning with some results for predictive guidance and then continuing to the velocity to be gained based guidance. In section IV simulations are presented to illustrate the results and in section V some conclusions are offered.

II.

Trajectory Model

Accurate prediction of trajectories for spinning projectiles is normally based on nonlinear models of the rigid body dynamics with six degrees of freedom (6DOF) and twelve state variables. For guidance purposes, however, details of the orientation dynamics are less important and an accurate description of the motion of the center of mass (CoM) is sufficient. Therefore, versions of the modified point mass model26 (MPM) is often used in guidance contexts. The MPM model has four degrees of freedom and seven state variables, or

3 of 24 American Institute of Aeronautics and Astronautics

if the decay in spin rate is predicted as a function of e.g. downrange, three degrees of freedom and six states. A model of the latter type, possibly with additional empirical corrections, can often yield predictions that are accurate enough for precision guidance and at the same time has a low complexity compared to the full 6DOF model. The guidance equations developed here are based on such a generic form of the equations of motion with three degrees of freedom and six state variables. II.A.

Firing scenario

Downloaded by John Robinson on January 26, 2015 | http://arc.aiaa.org | DOI: 10.2514/6.2015-0087

We shall consider firing of projectiles over a flat Earth. In the developments that follows the frame of reference can be arbitrary (nonrotating), but it is convenient to think of it as the fire control frame. The fire control frame is a Cartesian right handed frame with the direction given by the nominal (i.e. without aiming errors) orientation of the gun barrel; the 𝑥-direction is pointing downrange, the 𝑦-direction pointing upwards and the 𝑧-direction pointing crossrange to the right. We assume that the muzzle of the gun is located at a point with coordinates 𝑥𝐺 ∈ R3 and the target (desired terminal point) is located at a (fixed) point given by 𝑥𝑇 ∈ R3 (neither 𝑥𝐺 nor 𝑥𝑇 need represent a point on the Earth’s surface). We further assume that the target location relative to the gun is such that a firing solution exists, which we define next. II.B.

Firing solution

Much of what follows centers around the dependence on initial values for solutions to ordinary differential equations (ODEs) and for this reason we shall use the mathematical formalism of flows to describe the projectile trajectories.25 II.B.1.

Firing flow

Given the position of the gun 𝑥𝐺 ∈ R3 , a firing flow Φ𝑥𝐺 is a 𝐶 1 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 a family of solutions to the differential equations for the motion of the CoM of the projectile when the type of simplified models discussed above are used. Given also the position of the target 𝑥𝑇 ∈ R3 , a firing flow Φ𝑥𝐺 is a firing solution flow ifc Φ𝑥𝐺 ((𝑥𝐺 , 𝑣 𝐺 ), 𝑡𝑓 ) = (𝑥𝑇 , 𝑣 𝑓 )

(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. (Thus, a nominal trajectory is the trajectory obtained from a firing solution.) 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 (𝑥𝐺 , 𝑣 𝐺 ) ∈ ℳ (at time 𝑡 = 0) passes through the target location 𝑥𝑇 (at time 𝑡 = 𝑡𝑓 ). Since the flow Φ𝑥𝐺 is continuously differentiable, the position-velocity pair (𝑥1 , 𝑣 1 ) for the projectile at some later time 𝑡1 > 0 is a continuously differentiable function27 of the initial condition values (𝑥0 , 𝑣 0 ). The firing flow formalism employed 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. 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. c Note

that 𝑣 𝐺 and 𝑡𝑓 , 𝑣 𝑓 in (1) are in general not unique; compare the high and low elevation solutions for non maximal range fire.

4 of 24 American Institute of Aeronautics and Astronautics

Downloaded by John Robinson on January 26, 2015 | http://arc.aiaa.org | DOI: 10.2514/6.2015-0087

Figure 1: A firing flow which is also a firing solution flow since the trajectory starting at (𝑥𝐺 , 𝑣 𝐺 ) passes through the point 𝑥𝑇 at some later time. Another trajectory with a different initial velocity is also indicated which does not pass through 𝑥𝑇 and is thus not a firing solution.

II.C.

State space representation and controlled dynamics

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 , represents the aerodynamic acceleration at free flight trim (in some atmospheric conditions), i.e. yaw of repose (with control surfaces in their nominal setting, corresponding to 𝑢 = 0), and possibly Coriolis acceleration.d The control actionse are represented in open loop by the acceleration term 𝑢 : ℐ → R3 , which is a 𝐶 1 function of time. (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) represents the net effect of the actuators. For a fuze with nose mounted canards this translates to the acceleration resulting from the altered trim state (away from yaw of repose) due to deflection of the canards (not merely the effect from canard control forces). Since 𝑎 is nonlinearly dependent on the position-velocity pair (𝑥, 𝑣) the system of differential equations in (2) is nonlinear but it has (at least locally) a unique solution for any initial condition (𝑥0 , 𝑣 0 ) ∈ ℳ because 𝑎 and 𝑢 are 𝐶 1 . We assume that the model is formulated so that the solution is defined for all 𝑡 ∈ ℐ. The family of solutions obtained by varying the initial conditions (for a given 𝑢) defines a flow Φ : ℳ × ℐ → ℳ. As remarked in Ref. 25, the model (2) is quite general for guidance since it can it can incorporate empirical corrections computed with a more accurate method, e.g. a full a 6DOF model. 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 d The controlled flow is time dependent but we shall assume that this has formally been dealt with by enlarging ℳ with an extra state in the standard fashion (without explicitly indicating this in the notation). Propulsion could also be included. However, since guidance is normally applied after propulsion termination, it need not be explicitly taken into account here. e 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.

5 of 24 American Institute of Aeronautics and Astronautics

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 in the nominal time interval ℐ which is considered as starting point for a certain part of the discussion.f II.C.2.

Representation as linear perturbation of the vacuum trajectory

Downloaded by John Robinson on January 26, 2015 | http://arc.aiaa.org | DOI: 10.2514/6.2015-0087

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) ¯ which can (That is, we first solve (2) as it stands and then use the solution to define a time varying quantity 𝑎 be used instead of 𝑎 + 𝑢 when solving (2); the result (for fixed initial conditions) of the two approaches will ¯ is (unlike 𝑎) only directly dependent on time, on the initial condition (𝑥0 , 𝑣 0 ) be the same.) The function 𝑎 ¯ is and on the control function 𝑢( · ) (both in the open and closed loop cases). Furthermore, the function 𝑎 (by 𝐶 1 dependence of solutions to ODEs with respect to initial conditions, cf. e.g. Thm. 2.10 in Ref. 27) 𝐶 1 in the initial condition (𝑥0 , 𝑣 0 ) and (by continuous dependence of solutions to ODEs with respect to the vector field, cf. e.g. Thm. 3.4 in Ref. 28) continuous in sup-norm in the control 𝑢. 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 [︃ 𝐴=

03×3 03×3

𝐼 3×3 03×3

]︃

[︃ ,

𝐵=

03×3 𝐼 3×3

]︃

[︃ ,

exp(𝜉𝐴) = 𝐼 6×6 + 𝜉𝐴 =

𝐼 3×3 03×3

𝜉𝐼 3×3 𝐼 3×3

]︃ , 𝜉 ∈ R.

(5)

The formula (4) is a special case of a relation for perturbed nonlinear systems due to Alekseev.29, 30 ¯ obtained for 𝑢 = 0 will be denoted 𝑎 ¯ ((𝑥0 , 𝑣 0 ), 𝑡), i.e. we set The acceleration 𝑎 ¯ ((𝑥0 , 𝑣 0 ), 𝑡) = 𝑎 ¯ ((𝑥0 , 𝑣 0 ), 𝑡; 0). 𝑎

(6)

For zero air density (and nonrotating Earth) the term 𝑎 in (2) is identically zero without propulsion 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 ), 𝑡) ⃒⃒ ∆𝑥 + ∆𝑣 + ∆𝑡 + 𝑜 ‖((∆𝑥, ∆𝑣), ∆𝑡)‖ . (7) + 𝑥=𝑥 𝑣=𝑣 𝑠=𝑡 0 0 𝜕𝑥 𝜕𝑣 𝜕𝑠 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 given by the three perturbation terms on the right in (7). The guidance principle we propose below is based on this simple propertyg applied to (4). f A time shift of the flow means that initial values must be reinterpreted but this can can easily be done by mapping them back to original initial values of the unshifted flow using the group property of the flow (cf. (25) below). g The standard method to obtain perturbation formulas is to introduce the sensitivity equation, which is a linear time varying differential equation for the dynamics linearized around a nominal solution trajectory, cf. e.g. sec. 3.3 in Ref. 28. The approach used here is a very direct way to obtain estimates of the sensitivity matrix without actually solving any differential equation.

6 of 24 American Institute of Aeronautics and Astronautics

II.D.1.

Coriolis acceleration.

The Coriolis contribution to 𝑎 is normally very small and will essentially vanish in the perturbation expressions for the fire solution flow. Hence, we shall neglect the Coriolis contribution and 𝑎 will thus subsequently represent only aerodynamic acceleration. 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 thath [︃ ]︃ ∫︀ 𝑡 0 ,𝑣 0 ),𝑠) 𝜕Φ((𝑥0 , 𝑣 0 ), 𝑡) 𝐼 3×3 + 0 (𝑡 − 𝑠) 𝜕 𝑎¯ ((𝑥𝜕𝑥 𝑑𝑠 ∫︀ 𝑡 𝜕 𝑎¯ ((𝑥0 ,𝑣0 ),𝑠) = (8) 𝜕𝑥 𝑑𝑠 𝜕𝑥 0 and

Downloaded by John Robinson on January 26, 2015 | http://arc.aiaa.org | DOI: 10.2514/6.2015-0087

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

[︃

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

(9)

As a rough approximation (but in practice surprisingly useful25 ) one can use the vacuum trajectory version of the expressions in (8) and (9), which is obtained by setting the integrals in (8) and (9) to zero. For more accurate approximations, the contribution from the integrals needs to be taken into account. Since the integrands in the integrals in (8) and (9) are defined in terms of the nominal trajectory, the integrals can be computed offline and stored in tabulated form in the guidance computer. For the effect of a perturbation in the terminal time we get similarly in the uncontrolled case that [︃ ]︃ 𝜕Φ((𝑥0 , 𝑣 0 ), 𝑡) [Φ((𝑥0 , 𝑣 0 ), 𝑡)]4:6 = . (10) 𝜕𝑡 ¯ ((𝑥0 , 𝑣 0 ), 𝑡) 𝑔+𝑎 II.E.

Approximations

We shall now develop useful approximate expressions for the terms on the right hand sides of (8) and (9) as well as for the effects of (certain types of) 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 ,

(11)

and accordingly set 𝑉 (𝑡) = ‖𝑣(𝑡)‖ (so that in particular 𝑉0 = ‖𝑣 0 ‖). II.E.1.

Subsonic drag-only approximation

The drag-only approximation for 𝑎 in (2) is defined by setting 𝑎((𝑥, 𝑣), 𝑡) = 𝑎𝑑 ((𝑥, 𝑣)) where 𝜌(𝑥)𝑆𝐶𝐷 𝑉 𝑣. (12) 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 (12) is often expressed in terms of the ballistic coefficient 𝛽𝑑 = 𝑚/(𝑆𝐶𝐷 ), here defined in units of mass per area. The relation (12) gives 𝑎𝑑 ((𝑥, 𝑣)) = −

)︀𝑇 𝑆𝐶𝐷 𝑉 (︀ 𝜕𝑎𝑑 ((𝑥, 𝑣)) =− 𝑣 ∇𝜌(𝑥) . (13) 𝜕𝑥 2𝑚 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.i With this approximation we have )︀ 𝜌(𝑥)𝑆𝐶𝐷 𝑉 (︀ 1 𝜕𝑎𝑑 ((𝑥, 𝑣)) =− 𝐼 3×3 + 2 𝑣𝑣 𝑇 . (14) 𝜕𝑣 2𝑚 𝑉 h Differentiation i This

under the integral signs is easy to justify here and in the following. is reasonable for subsonic flight, cf. e.g.,26 for small aerodynamic angles.

7 of 24 American Institute of Aeronautics and Astronautics

II.E.2.

Functional approximations for the perturbation terms

In Ref. 25 it is shown that for a standard firing scenario for a 155mm artillery projectile and times up to 25s the upper part of the matrix on the right in (8) can be approximated as )︀𝑇 𝜕[Φ((𝑥0 , 𝑣 0 ), 𝑡)]1:3 𝑡2 𝑆𝐶𝐷 𝑉0 (︀ = 𝐼 3×3 − 𝑣 0 ∇𝜌(𝑥0 ) . 𝜕𝑥 2 2𝑚

(15)

For the upper part of the matrix on the right in (9) we have similarly the approximation )︀ 𝜕[Φ((𝑥0 , 𝑣 0 ), 𝑡)]1:3 𝑡2 𝜌(𝑥0 )𝑆𝐶𝐷 𝑉0 (︀ 1 = 𝑡𝐼 3×3 − 𝐼 3×3 + 2 𝑣 0 𝑣 𝑇0 . 𝜕𝑣 2 2𝑚 𝑉0 II.F.

(16)

Effect of controls

Downloaded by John Robinson on January 26, 2015 | http://arc.aiaa.org | DOI: 10.2514/6.2015-0087

The solution Φ((𝑥0 , 𝑣 0 ), 𝑡; 𝑢) in (4) is nonlinear in the control 𝑢. However, for small 𝑢 the influence can be approximated as linear and we shall develop one such approximation here, for a certain class of controls. II.F.1.

Form of the control

Depending on the type of guidance considered we make different assumptions about the control 𝑢. In the developments below on guidance based on velocity to be gained techniques we assume that the control 𝑢 is a ¯ will be continuously general 𝐶 1 function. From the discussion above we know that the values of both 𝑎 and 𝑎 dependent on 𝑢 (sup-norm). In the developments below of predictive guidance methods we assume that the control function 𝑢 is a pulse function and of the form 𝑢(𝑡) = 𝐻(𝑡)𝑢𝜃 , (17) where 𝐻 : ℐ → R3×3 is a matrix valued 𝐶 1 function and 𝑢𝜃 ∈ R3 is a parameter vector. By 𝐶 1 dependence ¯ in (3) are continuously differentiable with on parameters for solutions to ODEs it follows that 𝑎 in (2) and 𝑎 respect to 𝑢𝜃 . We shall need to make also some further assumptions about the control 𝑢 in (17). Assumptions about the pulse function. We make the following three assumptions about 𝐻 in (17). (i): The matrix 𝐻(𝑡) is nonsingular, for all 𝑡 ∈ ℐ. With this assumption, there is locally a bijective relation between the values 𝑢(𝑡) of the control and the parameter vector 𝑢𝜃 . (ii): The integral ∫︁ 𝑡 (𝑡 − 𝑠)𝐻(𝑠) 𝑑𝑠, 0

is nonsingular for 𝑡 ∈ ℐ, except possibly a finite number of points. This assumption will guarantee that there is a control (guidance command) that solves our guidance problem below, and that this control is unique.j (iii): The matrix 𝐻 has been normalized such that ‖𝐻(𝑡)𝜂‖ = 1. ‖𝜂‖ 𝑡∈ℐ,𝜂∈R3 ∖{0} sup

(18)

This third assumption is introduced to ensure that the components of 𝑢𝜃 really represent amplitudes of the control signals so that some estimates below behave the same way as they would in case 𝑢(𝑡) = 𝑢𝜃 , i.e. if 𝑢 were constant. In real applications the function 𝐻 is generally redefined in each guidance update, e.g. by time shifting/stretching a basic pulse function. Therefore, in the formulas below where the time 𝑡 = 0 can refer to any starting point in the nominal time interval ℐ the proper shift of 𝐻 has to be inferred from the context. j It

is easy to see that assumption (i) is not sufficient to guarantee (ii) and that both assumptions can be expressed as the ∫︀ assumption that 𝑡 ↦→ 0𝑡 ℎ(𝑡 − 𝑠)𝐻(𝑠) 𝑑𝑠 is nonsingular in general, for ℎ(𝑡 − 𝑠) = 1 and ℎ(𝑡 − 𝑠) = 𝑡 − 𝑠, respectively.

8 of 24 American Institute of Aeronautics and Astronautics

Basic estimate of the influence of control. When the control 𝑢 is of the form (17) we have by 𝐶 1 -smoothness of 𝑎, Φ thatk (︀ )︀ (︀ )︀ (︀ )︀ ¯ (𝑥0 , 𝑣 0 ), 𝑡; 𝑢) 𝜕𝑎 𝜕𝑎 Φ((𝑥0 , 𝑣 0 ), 𝑡; 𝑢), 𝑡 𝜕𝑎 Φ((𝑥0 , 𝑣 0 ), 𝑡; 𝑢), 𝑡 𝜕Φ((𝑥0 , 𝑣 0 ), 𝑡; 𝑢) = + 𝐻(𝑡) = + 𝐻(𝑡) 𝜕𝑢𝜃 𝜕𝑢𝜃 𝜕(𝑥, 𝑣) 𝜕𝑢𝜃 (︀ )︀ (︀ )︀ 𝜕𝑎 Φ((𝑥0 , 𝑣 0 ), 𝑡; 𝑢), 𝑡 𝜕[Φ((𝑥0 , 𝑣 0 ), 𝑡; 𝑢)]1:3 𝜕𝑎 Φ((𝑥0 , 𝑣 0 ), 𝑡; 𝑢), 𝑡 𝜕[Φ((𝑥0 , 𝑣 0 ), 𝑡; 𝑢)]4:6 = + + 𝐻(𝑡). 𝜕𝑥 𝜕𝑢𝜃 𝜕𝑣 𝜕𝑢𝜃 (19) For a vacuum trajectory (𝑎 = 0) we have

Vacuum trajectory approximation.

¯ ((𝑥0 , 𝑣 0 ), 𝑡; 0)) 𝜕𝑎 = 𝐻(𝑡). 𝜕𝑢𝜃

(20)

Downloaded by John Robinson on January 26, 2015 | http://arc.aiaa.org | DOI: 10.2514/6.2015-0087

Drag-only approximation. To investigate more accurate approximations for the right hand side of (19) based on the drag-only approximation (12) we recall (4) and (5). From these relations it follows that 𝜕Φ((𝑥0 , 𝑣 0 ), 𝑡; 0) = 𝜕𝑢𝜃

∫︁ 0

𝑡

(︀ )︀ 𝜕 𝑎 ¯ ((𝑥0 , 𝑣 0 ), 𝑠; 0)) exp (𝑡 − 𝑠)𝐴 𝐵 𝑑𝑠 = 𝜕𝑢𝜃 ∫︁ 𝑡 [︃ 0

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

]︃

¯ ((𝑥0 , 𝑣 0 ), 𝑠; 0)) 𝜕𝑎 𝑑𝑠. 𝜕𝑢𝜃

Combining this with (19) and (6) yields the functional equation ¯ ((𝑥0 , 𝑣 0 ), 𝑡; 0)) ¯ ((𝑥0 , 𝑣 0 ), 𝑡; 0) 𝜕𝑎 𝜕𝑎 = 𝜕𝑢𝜃 𝜕𝑥

∫︁

𝑡

¯ ((𝑥0 , 𝑣 0 ), 𝑠; 0)) 𝜕𝑎 𝑑𝑠 𝜕𝑢𝜃 ∫︁ ¯ ((𝑥0 , 𝑣 0 ), 𝑠; 0)) ¯ ((𝑥0 , 𝑣 0 ), 𝑡; 0) 𝑡 𝜕 𝑎 𝜕𝑎 + 𝑑𝑠 + 𝐻(𝑡). (21) 𝜕𝑣 𝜕𝑢𝜃 0

(𝑡 − 𝑠) 0

¯ /𝜕𝑢𝜃 where the (integral) operator in question is contractive for This is an affine functional equation for 𝜕 𝑎 sufficiently small 𝑡 and then the equation has a unique solution. In Ref. 25 it is shown that for a standard 155mm artillery projectile firing scenario and times 𝑡 up to (at least) about 10s it is a good approximation to use the vacuum trajectory approximation in (20) as the solution to (21) when 𝐻(𝑡) = 𝐼 3×3 . By a slight extension of that argument it can be shown that the vacuum trajectory solution to (21) is a good approximation also for a 𝐻(𝑡) as in (17), provided (18) holds. Therefore, we shall use the vacuum trajectory ¯ in (6) of control when the control 𝑢 is of approximation in (20) as the basic estimate for the influence on 𝑎 the form (17). II.F.2.

Perturbation expression for control

With the simple approximation (20), the effect of a small variation (around zero) of control parameters 𝑢𝜃 in a control 𝑢 as in (17) is easily calculated from the variation-of-parameters formula (4) and we obtain ]︃ [︃ ∫︀ 𝑡 ]︃ [︃ ∫︀ ¯ ((𝑥0 ,𝑣 0 ),𝑠;0) 𝜕𝑎 𝑡 (𝑡 − 𝑠) 𝑑𝑠 𝜕Φ((𝑥0 , 𝑣 0 ), 𝑡; 0) (𝑡 − 𝑠)𝐻(𝑠) 𝑑𝑠 𝜕𝑢𝜃 0 ∫︀ 0 ∫︀ = = . 𝑡 𝑡 𝜕𝑎 ¯ ((𝑥0 ,𝑣 0 ),𝑠;0) 𝜕𝑢𝜃 𝑑𝑠 𝐻(𝑠) 𝑑𝑠 𝜕𝑢𝜃 0 0 A first order approximation of the effect in Φ in (4) of nonzero control parameters 𝑢𝜃 is thus given by 𝜕Φ((𝑥0 , 𝑣 0 ), 𝑡; 0) 𝑢𝜃 + 𝑜(‖𝑢𝜃 ‖) = 𝜕𝑢𝜃 ∫︁ 𝑡 ∫︁ 𝑡 )︀ (︀ Φ((𝑥0 , 𝑣 0 ), 𝑡) + (𝑡 − 𝑠)𝐻(𝑠)𝑢𝜃 𝑑𝑠, 𝐻(𝑠)𝑢𝜃 𝑑𝑠 + 𝑜(‖𝑢𝜃 ‖). (22)

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

0

0

This is the basic relation for the influence of control that we shall use for controls 𝑢 of the form (17). It is a generalization to pulse function control of the corresponding relation for constant control given in Ref. 25. k Smoothness

of the solution Φ with respect to 𝑢𝜃 follows from standard properties of flows (continuity with respect to initial values) since 𝑢𝜃 can be regarded as an extra state variable with zero time derivative.

9 of 24 American Institute of Aeronautics and Astronautics

III.

Guidance

We shall now use the perturbation expressions (7)–(10) (with the approximations in (15), (16)) to develop equations for guidance. The main result is a closed loop guidance scheme based on the concept of velocity to be gained. The velocity to be gained is calculated via the predicted zero effort terminal deviation and we start by deriving a predictive guidance scheme based on this quantity. III.A.

Predictive guidance

Let Φ𝑥𝐺 be a firing solution flow over a time interval ℐ, as in section II.B, with firing solution (𝑥𝐺 , 𝑥𝑇 , 𝑣 𝐺 ) and nominal trajectory (with control set to zero, 𝑢 = 0) 𝑡 ↦→ Φ𝑥𝐺 ((𝑥𝐺 , 𝑣 𝐺 ), 𝑡),

𝑡 ≥ 0,

(23)

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

(24)

Downloaded by John Robinson on January 26, 2015 | http://arc.aiaa.org | DOI: 10.2514/6.2015-0087

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

(25)

where 𝑡1 is any time point in the nominal time-of-flight interval [0, 𝑡𝑇 ]. 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 (24) and the group property (25), the deviation 𝜁 0 (𝑡𝐹 ) ∈ R6 between the terminal state Φ𝑥𝐺 ((𝑥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 ). (26) The deviation 𝜁 0 (𝑡𝐹 ) in (26) will be called the zero effort (terminal) deviation (at time 𝑡𝐹 ) (ZED). If we apply the first order perturbation approximation in (7) to 𝜁 0 (𝑡𝐹 ) we obtain, approximately (for small ∆𝑥, ∆𝑣, ∆𝑡) 𝜕Φ𝑥𝐺 (Φ𝑥𝐺 ((𝑥𝐺 , 𝑣 𝐺 ), 𝑡1 ), 𝑡𝑇 − 𝑡1 )∆𝑥 𝜕𝑥 𝜕Φ𝑥𝐺 𝜕Φ𝑥𝐺 + (Φ𝑥𝐺 ((𝑥𝐺 , 𝑣 𝐺 ), 𝑡1 ), 𝑡𝑇 − 𝑡1 )∆𝑣 + (Φ𝑥𝐺 ((𝑥𝐺 , 𝑣 𝐺 ), 𝑡1 ), 𝑡𝑇 − 𝑡1 )∆𝑡, 𝜕𝑣 𝜕𝑡

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

(27)

where ∆𝑥 = 𝑥1 − [Φ𝑥𝐺 ((𝑥𝐺 , 𝑣 𝐺 ), 𝑡1 )]1:3 , ∆𝑣 = 𝑣 1 − [Φ𝑥𝐺 ((𝑥𝐺 , 𝑣 𝐺 ), 𝑡1 )]4:6 , ∆𝑡 = 𝑡𝐹 − 𝑡𝑇 . III.A.2.

Terminal position deviation with control

It is straightforward to compute an approximation to the deviation in terminal state obtained when a (small) control acceleration 𝑢 of the form (17) is applied to the firing flow, using the relation (22) and the ZED expressions above. With notation as in section III.A.1, consider the controlled terminal deviation (CTD) 𝜁 𝑢 (𝑡𝐹 ) ∈ R3 at time 𝑡𝐹 defined by 𝜁 𝑢 (𝑡𝐹 ) = Φ𝑥𝐺 ((𝑥1 , 𝑣 1 ), 𝑡𝐹 − 𝑡1 ; 𝑢) − (𝑥𝑇 , 𝑣 𝑇 )

10 of 24 American Institute of Aeronautics and Astronautics

(28)

Downloaded by John Robinson on January 26, 2015 | http://arc.aiaa.org | DOI: 10.2514/6.2015-0087

Figure 2: Closest point of approach (CPA). The actual flight path is solid and the nominal is dashed. Evaluating the error along the actual flight path but at the time of nominal impact may lead to inaccurate results. For example, if the actual flight path is nearly a “delayed” or “advanced” version of the nominal, an apparent terminal error would be present which in reality might be small and not require correction.

obtained when a 𝐶 1 control acceleration 𝑢 is applied to the perturbed trajectory. From (22) we then have that the deviation in terminal state with control of the form (17) is approximately (for small 𝑢) 𝜁 𝑢 (𝑡𝐹 ) = Φ𝑥𝐺 ((𝑥1 , 𝑣 1 ), 𝑡𝐹 − 𝑡1 ; 𝑢) − (𝑥𝑇 , 𝑣 𝑇 ) = ∫︁ ∫︁ 𝑡𝐹 −𝑡1 (︀ 𝑡𝐹 −𝑡1 )︀ Φ𝑥𝐺 ((𝑥1 , 𝑣 1 ), 𝑡𝐹 − 𝑡1 ) − (𝑥𝑇 , 𝑣 𝑇 ) + (𝑡𝐹 − 𝑡1 − 𝑠)𝐻(𝑠)𝑢𝜃 𝑑𝑠, 𝐻(𝑠)𝑢𝜃 𝑑𝑠 = 0 0 ∫︁ ∫︁ 𝑡𝐹 −𝑡1 (︀ 𝑡𝐹 −𝑡1 )︀ 𝜁 0 (𝑡𝐹 ) + (𝑡𝐹 − 𝑡1 − 𝑠)𝐻(𝑠)𝑢𝜃 𝑑𝑠, 𝐻(𝑠)𝑢𝜃 𝑑𝑠 , (29) 0

0

where 𝜁 0 (𝑡𝐹 ) is the ZED in (26). Since the state vector (𝑥, 𝑣) is six dimensional but the control 𝑢 is only three dimensional (when 𝐻 in (17) is fixed) there is in general no solution for the problem of driving the entire state deviation 𝜁 𝑢 (𝑡𝐹 ) in (29) to zero during [𝑡𝐹 − 𝑡1 , 𝑡𝑇 ] using such a control 𝑢 (regardless of 𝑡𝐹 ). However, the position part of the state can be nulled and for this we can use (29) to obtain approximately (for small 𝑢) that ∫︁ 𝑡𝐹 −𝑡1 [𝜁 𝑢 (𝑡𝐹 )]1:3 = 𝜁 0,𝑥 (𝑡𝐹 ) + (𝑡𝐹 − 𝑡1 − 𝑠)𝐻(𝑠)𝑢𝜃 𝑑𝑠, (30) 0 3

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

(31)

Two dimensional error manifolds. The ZEPD 𝜁 0,𝑥 (𝑡𝐹 ) in (31) takes values in R3 but in our man application (2D course correcting fuzes) the control takes values in a two-dimensional manifold (lateral and longitudinal acceleration). Therefore, it is of interest to reduce the dimension of the manifold on which the error is defined so that there can be one-to-one correspondence between error and control variables. There are several ways to do this but a straightforward way is to simply project the error onto the plane perpendicular to the predicted terminal velocity, cf. figure 2, which can be motivated as follows. When the projectile is close to the target, and flies essentially along a straight path, the terminal error in the direction along the flight path can largely be ignored since it essentially corresponds to a change in impact time. 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 ‖𝜁 0,𝑥 (𝑡𝐹 )‖ ≥ ‖𝜁 0,𝑥 (𝑡𝐶𝑃 𝐴 )‖,

𝑡𝐹 ∈ ℐ,

(32)

where 𝑡𝐶𝑃 𝐴 ∈ ℐ is the time such that the orthogonality conditionl 𝜁 0,𝑥 (𝑡𝐶𝑃 𝐴 )𝑇 [Φ𝑥𝐺 ((𝑥1 , 𝑣 1 ), 𝑡𝐶𝑃 𝐴 − 𝑡1 )]4:6 = 0 l We

assume that, in some neighborhood of 𝑡𝑇 , there is precisely one time 𝑡𝐶𝑃 𝐴 such that (33) holds.

11 of 24 American Institute of Aeronautics and Astronautics

(33)

is satisfied, see fig. 2. Using the vacuum trajectory approximation (cf. (8), (9)) the position [Φ𝑥𝐺 ((𝑥1 , 𝑣 1 ), 𝑡𝐹 − 𝑡1 )]1:3 (at an arbitrary terminal time 𝑡𝐹 ) can be estimated as [Φ𝑥𝐺 ((𝑥1 , 𝑣 1 ), 𝑡𝐹 − 𝑡1 )]1:3 = [Φ𝑥𝐺 ((𝑥1 , 𝑣 1 ), 𝑡𝑇 − 𝑡1 )]1:3 + (𝑡𝐹 − 𝑡𝑇 )[Φ𝑥𝐺 ((𝑥1 , 𝑣 1 ), 𝑡𝑇 − 𝑡1 )]4:6 , and an associated estimate for 𝜁 0,𝑥 (𝑡𝐶𝑃 𝐴 ) follows from (31) as 𝜁 0,𝑥 (𝑡𝐶𝑃 𝐴 ) = 𝜁 0,𝑥 (𝑡𝑇 ) + (𝑡𝐶𝑃 𝐴 − 𝑡𝑇 )[Φ𝑥𝐺 ((𝑥1 , 𝑣 1 ), 𝑡𝑇 − 𝑡1 )]4:6 .

(34)

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

Downloaded by John Robinson on January 26, 2015 | http://arc.aiaa.org | DOI: 10.2514/6.2015-0087

𝑡𝐶𝑃 𝐴 − 𝑡𝑇 = − and

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

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

(35)

where 𝑃 [𝑣𝑇 ]⊥ denotes projection onto the subspace orthogonal to 𝑣 𝑇 . Driving the terminal position deviation to zero. If we take 𝑡𝐹 = 𝑡𝐶𝑃 𝐴 in (30) we obtain a simple approximate relation (for small 𝑢) that can be used to chose a control acceleration 𝑢 of the form (17) to drive the predicted position part [𝜁 𝑢 (𝑡𝐶𝑃 𝐴 )]1:3 of the controlled terminal deviation to (approximately) zero ˆ 𝜃 where 𝑢 ˆ 𝜃 is given by during [𝑡1 , 𝑡𝐶𝑃 𝐴 − 𝑡1 ]. The solution is to set 𝑢𝜃 = 𝑢 ˆ𝜃 = − 𝑢

(︁ ∫︁

𝑡𝐶𝑃 𝐴 −𝑡1

)︁−1 (𝑡𝐶𝑃 𝐴 − 𝑡1 − 𝑠)𝐻(𝑠) 𝑑𝑠 𝜁 0,𝑥 (𝑡𝐶𝑃 𝐴 )

(36)

0

and 𝜁 0,𝑥 (𝑡𝐶𝑃 𝐴 ) is the (uncontrolled) CPA in (32), which can be approximated as in (35). The inverse on the right (36) is well-defined, at least locally around 𝑢𝜃 = 0, due to assumption (ii) on 𝐻 in sec. II.F.1. III.A.3.

Predictive perturbation guidance algorithm

A predictive guidance algorithm can now be given as follows. Before flight, compute a nominal trajectory (23) using any kind of high fidelity method. Then, at each guidance update instant during flight, do: (1) Compute an approximation to the ZEPD 𝜁 0,𝑥 (𝑡𝐹 ) in (31) for 𝑡𝐹 = 𝑡𝑇 (the nominal terminal time). The quantity [Φ𝑥𝐺 ((𝑥1 , 𝑣 1 ), 𝑡𝑇 − 𝑡1 )]1:3 is then approximated by setting 𝑥1 = 𝑥0 + ∆𝑥, 𝑣 1 = 𝑣 0 + ∆𝑣, 𝑡 = 𝑡𝑇 − 𝑡1 and considering the three perturbation terms on the right in (7). The Jacobians in the perturbation terms are obtained from (8), (9) where the approximations in section II.E.2 are used. (2) Compute approximations of 𝑡𝐶𝑃 𝐴 − 𝑡𝑇 and 𝜁 0,𝑥 (𝑡𝐶𝑃 𝐴 ) in (32) using the formula in (35). ˆ 𝜃 given by (36). (3) Compute the commanded acceleration 𝑢(𝑡) as in (17) with 𝑢𝜃 = 𝑢 This algorithm generalizes the one given in Ref. 25. Both algorithms are inherently open loop schemes but can be used in closed loop by performing the computations repeatedly. The advantage with having pulse functions as in (17) for control is that the prediction and actuation intervals can be separated. III.B.

Velocity to be gained guidance

We now continue and develop a closed loop guidance scheme based on the concept of velocity to be gained. The scheme is obtained by applying feedback to the open loop dynamics for the velocity to be gained, which we derive first.

12 of 24 American Institute of Aeronautics and Astronautics

III.B.1.

Required velocity

Consider a firing solution flow Φ𝑥𝐺 with nominal trajectory as in (23) and terminal condition (24). Let 𝑡0 ∈ (0, 𝑡𝑇 ) be arbitrary, assume that the control is zero 𝑢 = 0 and put (𝑥0 , 𝑣 0 ) = Φ𝑥𝐺 ((𝑥𝐺 , 𝑣 𝐺 ), 𝑡0 ),

(37)

i.e. (𝑥0 , 𝑣 0 ) denotes a point along a nominal trajectory (a trajectory leading to the target at 𝑥𝑇 ). Assume that there exist neighborhoods 𝑈 ⊆ R3 , 𝑉 ⊆ R3 with 𝑥0 ∈ 𝑈, 𝑣 0 ∈ 𝑉 such that the equationm

Downloaded by John Robinson on January 26, 2015 | http://arc.aiaa.org | DOI: 10.2514/6.2015-0087

[Φ𝑥𝐺 ((𝑥, 𝑣), 𝑡𝑇 − 𝑡0 )]1:3 = 𝑥𝑇

(38)

has a unique solution 𝑣 ∈ 𝑉 for any 𝑥 ∈ 𝑈 . If we let 𝑥 be the position at time 𝑡0 on a perturbed trajectory in the flow Φ𝑥𝐺 it is clear that the solution 𝑣 can be interpreted as the velocity that is needed at (𝑥, 𝑡0 ) in order to reach the nominal terminal position 𝑥𝑇 at time 𝑡𝑇 . Given a time 𝑡0 ∈ (0, 𝑡𝑇 ), the map 𝑣 𝑟 : 𝑈 → 𝑉 defined by 𝑥 ↦→ 𝑣 via the solution to (38) will be called the required velocity and we note that in particular 𝑣 𝑟 (𝑥0 , 𝑡0 ) = 𝑣 0 , i.e. all velocities along the nominal trajectory are also required velocities for the corresponding positions. A sufficient condition for the required velocity 𝑣 𝑟 to be well defined at time 𝑡0 in a neighborhood of the position 𝑥0 on the nominal trajectory (37) is that the Jacobian 𝜕[Φ𝑥𝐺 ((𝑥0 , 𝑣 0 ), 𝑡𝑇 − 𝑡0 )]1:3 𝜕𝑣

(39)

has full rank since the implicit function theorem then guarantees that neighborhoods 𝑈, 𝑉 as in (38) exist with 𝑥0 ∈ 𝑈, 𝑣 0 ∈ 𝑉 . Glancing back at (9) we see that it is reasonable to assume that a similar property holds more globally in time since we have the approximation in (16). Assumption about the position-velocity sensitivity. We shall assume that there exist two (time invariant) neighborhoods 𝑈0 , 𝑉0 ⊆ R3 containing 0 such that 𝜕[Φ𝑥𝐺 ((𝑥1 , 𝑣 1 ), 𝑡𝑇 − 𝑡0 )]1:3 , 𝜕𝑣

(40)

has full rank whenever (𝑥1 − 𝑥0 , 𝑣 1 − 𝑣 0 ) ∈ 𝑈0 × 𝑉0 and 𝑡0 ∈ (0, 𝑡𝑇 ), where (𝑥0 , 𝑣 0 ) is given by (37). When this assumptionn is fulfilled there exists “box tube” of position-velocity pairs around the nominal trajectory (23) where required velocity 𝑣 𝑟 is well defined.o Since the flow Φ𝑥𝐺 is 𝐶 1 in all arguments it follows by the implicit function theorem (smooth version) that 𝑥 ↦→ 𝑣 𝑟 (𝑥, 𝑡0 ) is 𝐶 1 too on 𝑈0 + {𝑥0 }. Required velocity vector field. Let (𝑥0 , 𝑡0 ) be a space-time point as in (37). The assumption relating to (40) guarantees that if 𝑥1 ∈ 𝑈0 + {𝑥0 } then the required velocity 𝑣 𝑟 is well defined and 𝐶 1 in a space-time neighborhood of (𝑥1 , 𝑡0 ). The required velocity 𝑣 𝑟 therefore represents a time dependent smooth vector field on the time varying manifold 𝑈0 + {𝑥0 } ⊆ R3 . The integral curves to this vector field are given by the solutions to ˙ 𝑥(𝑡) = 𝑣 𝑟 (𝑥(𝑡), 𝑡). (41) This differential equation has well defined local solutions (parametrized by initial conditions) but we have also the following global property which is key to the subsequent results. Proposition III.1. The solutions to (41) are given by (𝑥(𝑡), 𝑣(𝑡)) = Φ𝑥𝐺 ((𝑥1 , 𝑣 𝑟 (𝑥1 , 𝑡0 )), 𝑡 − 𝑡0 ),

𝑡 ∈ [𝑡0 , 𝑡𝑇 ),

(42)

where (𝑥1 , 𝑡0 ) is any space-time point such that 𝑥1 ∈ 𝑈0 + {𝑥0 }. m Like 𝑣 𝐺 in the firing solution, we can only expect a required velocity to be locally unique. In the case of powered flight the solution is in generally not even locally unique and additional constraints must be enforced to make it unique, cf. e.g. Ref. 31. n With the size estimates in Ref. 25 applied to (16) the assumption in (40) is fulfilled for 𝑡 − 𝑡 ≤ 80s (approximately). 0 𝑇 o Near the right end point of the interval [0, 𝑡 ] the Jacobian (39) is of the form (𝑡 − 𝑡 )𝐼 + 𝒪((𝑡 − 𝑡 )2 ). Hence, the 0 0 𝑇 𝑇 𝑇 Jacobian is guaranteed to have full rank when the time-to-go 𝑡𝑇 − 𝑡0 is sufficiently small but it also vanishes with the rate 𝑡𝑇 − 𝑡0 . Consequently, 𝑣 𝑟 can be expected to grow inversely with 𝑡𝑇 − 𝑡0 in 𝑈0 + {𝑥0 } when 𝑡𝑇 − 𝑡0 goes to 0.

13 of 24 American Institute of Aeronautics and Astronautics

Proof. Assume that 𝑥1 ∈ 𝑈0 + {𝑥0 } for some 𝑡0 < 𝑡𝑇 . Then, the definition of required velocity gives [Φ𝑥𝐺 ((𝑥1 , 𝑣 𝑟 (𝑥1 , 𝑡0 )), 𝑡𝑇 − 𝑡0 )]1:3 = 𝑥𝑇 .

(43)

If we define (𝑥(𝑡), 𝑣(𝑡)) as in (42) and apply (43) together with the flow group property (25) we have [Φ𝑥𝐺 ((𝑥(𝑡), 𝑣(𝑡)), 𝑡𝑇 − 𝑡)]1:3 = [Φ𝑥𝐺 ((𝑥1 , 𝑣 𝑟 (𝑥1 , 𝑡0 )), 𝑡𝑇 − 𝑡0 )]1:3 = 𝑥𝑇 ,

𝑡 ∈ [𝑡0 , 𝑡𝑇 ).

If we compare with the definition (38) we see that for each pair (𝑥(𝑡), 𝑣(𝑡)) as in (42) the velocity 𝑣(𝑡) equals the required velocity at the position 𝑥(𝑡) and since the flow Φ𝑥𝐺 represents solutions to (2) we have at the ˙ same time 𝑥(𝑡) = 𝑣(𝑡), i.e. (41) is satisfied. By uniqueness of solutions to (41) the result now follows. The result shows that if required velocity is attained 𝑣 = 𝑣 𝑟 at some time 𝑡0 < 𝑡𝑇 for a solution to (2) then the position variable for that solution will, if the control is held at zero 𝑢 = 0, continue and coast along Φ𝑥𝐺 towards 𝑥𝑇 during [𝑡0 , 𝑡𝑇 ). An immediate consequence of the result above is the following.

Downloaded by John Robinson on January 26, 2015 | http://arc.aiaa.org | DOI: 10.2514/6.2015-0087

Corollary III.2. Let (𝑥1 , 𝑡0 ) be a space-time point such that 𝑥1 ∈ 𝑈0 + {𝑥0 }. Then, the required velocity 𝑣 𝑟 satisfies (︀ )︀ 𝑑 𝜕𝑣 𝑟 (𝑥(𝑡), 𝑡) 𝜕𝑣 𝑟 (𝑥(𝑡), 𝑡) 𝑣 𝑟 (𝑥(𝑡), 𝑡) = + 𝑣 𝑟 (𝑥(𝑡), 𝑡) = 𝑔 + 𝑎 𝑥(𝑡), 𝑣 𝑟 (𝑥(𝑡), 𝑡) , 𝑑𝑡 𝜕𝑡 𝜕𝑥

𝑡 ∈ [𝑡0 , 𝑡𝑇 ),

(44)

where 𝑥(𝑡) is given by (42) (and satisfies (41)) with initial condition 𝑥(𝑡0 ) = 𝑥1 . The 𝑄-matrix. Again, let (𝑥0 , 𝑡0 ) be a space-time point as in (37), so that 𝑣 𝑟 is a 𝐶 1 vector field on 𝑈0 + {𝑥0 }. From the definition of 𝑣 𝑟 in (38) we have [Φ𝑥𝐺 ((𝑥, 𝑣 𝑟 (𝑥, 𝑡)), 𝑡𝑇 − 𝑡0 )]1:3 = 𝑥𝑇 ,

𝑥 ∈ 𝑈0 + {𝑥0 },

(45)

and it follows that 𝜕[Φ𝑥𝐺 ((𝑥, 𝑣 𝑟 (𝑥, 𝑡)), 𝑡𝑇 − 𝑡0 )]1:3 𝜕[Φ𝑥𝐺 ((𝑥, 𝑣 𝑟 (𝑥, 𝑡)), 𝑡𝑇 − 𝑡0 )]1:3 + 𝑄(𝑥, 𝑡0 ) = 0, 𝜕𝑥 𝜕𝑣

𝑥 ∈ 𝑈0 + {𝑥0 },

(46)

where

𝜕𝑣 𝑟 (𝑥, 𝑡0 ) . (47) 𝜕𝑥 The matrix in (47) is often called the 𝑄-matrix in space vehicle applications and will be used several times below. Since the the rightmost Jacobian in (46) is nonsingular by the assumption in (40) we can write 𝑄(𝑥, 𝑡0 ) =

(︀ 𝜕[Φ𝑥𝐺 ((𝑥, 𝑣 𝑟 (𝑥, 𝑡0 )), 𝑡𝑇 − 𝑡0 )]1:3 )︀−1 𝜕[Φ𝑥𝐺 ((𝑥, 𝑣 𝑟 (𝑥, 𝑡0 )), 𝑡𝑇 − 𝑡0 )]1:3 𝑄(𝑥, 𝑡0 ) = − . 𝜕𝑣 𝜕𝑥

(48)

It is of interest to obtain approximations for the right hand side in (48). By recalling the approximations in (15) and (16) we see that the dependence on both the position and velocity component in these expressions is weak (and is asymptotically zero as the time-to-go 𝑡𝑇 − 𝑡0 tends to 0). For 𝑥 ∈ 𝑈0 + {𝑥0 } we can therefore, in order to obtain an approximation, substitute (𝑥0 , 𝑣 0 ) for (𝑥, 𝑣 𝑟 (𝑥, 𝑡0 )) in (48). If we do this and use (15) and (16), and apply a geometric (Neumann) series operator expansion (validp for small 𝑡𝑇 − 𝑡0 > 0) we obtain the approximationq 𝑄(𝑥, 𝑡0 ) = −

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

𝑥 ∈ 𝑈0 + {𝑥0 }.

(49)

In the vacuum trajectory case 𝑎 = 0 (which here also implies constant gravity) we have (cf. (8), (9)) the well-known relation 𝑄(𝑥, 𝑡0 ) = 𝑄0 (𝑥, 𝑡0 ) where 𝑄0 (𝑥, 𝑡0 ) = −

1 𝐼 3×3 , 𝑡𝑇 − 𝑡0

p The

(50)

expansion exists when the assumption (40) can be expected to hold, see earlier footnote. this approximation, 𝑄 is dependent only on the nominal position 𝑥0 at 𝑡0 . In space vehicle applications it is the nonconstant gravity that gives rise to the 𝑥-dependence in 𝑄 and even though the dependence is weak, it is important for longer flight times. In our application, with constant gravity and moderately long flight times, the 𝑥-dependence in 𝑄 enters via aerodynamics (of the future flight) and is also weak. Therefore, it can be approximated by a dependence only on 𝑥0 . q With

14 of 24 American Institute of Aeronautics and Astronautics

Downloaded by John Robinson on January 26, 2015 | http://arc.aiaa.org | DOI: 10.2514/6.2015-0087

Figure 3: Velocity to be gained (VTG). The solid thick curve is a nominal trajectory leading to the target, the dashed is a perturbed trajectory not leading to the target and the thin solid is a trajectory leading to the target. The velocity to be gained 𝑣 𝑔 = 𝑣 𝑟 − 𝑣 is the difference between the required velocity 𝑣 𝑟 needed to reach the target without further controls (thin solid curve) and the actual (perturbed) velocity 𝑣.

and the 𝑄-matrix is then independent of the position. However, it is dependent on the time-to-go 𝑡𝑇 − 𝑡0 and becomes singular when 𝑡0 approaches 𝑡𝑇 . (This is inevitable regardless of the form of 𝑎 in (2).) So far we have defined the required velocity 𝑣 𝑟 and established some of its properties, including approximations of related quantities, but we have not given an explicit (easily computable) expression for 𝑣 𝑟 itself. This will be given in terms of the zero effort terminal position deviation, after we have introduced the velocity to be gained. We note also that everything above, including the required velocity, is defined for the case of zero control. The effect of nonzero control will be included in the definition of velocity to be gained. III.B.2.

Velocity to be gained

Let 𝑡 ↦→ (𝑥(𝑡), 𝑣(𝑡)) be any solution to (2), where the control 𝑢 is a 𝐶 1 function (not necessarily parametrized as in (17)). For each 𝑡0 ∈ (0, 𝑡𝑇 ) and 𝑥 ∈ 𝑈0 + {𝑥0 } define the velocity to be gained (VTG) 𝑣 𝑔 as 𝑣 𝑔 (𝑥, 𝑡0 ) = 𝑣 𝑟 (𝑥, 𝑡0 ) − 𝑣(𝑡0 ).

(51)

The velocity to be gained 𝑣 𝑔 can be interpreted as a time dependent vector field, defined in terms of (the velocity values 𝑣(𝑡0 ) of) one particular (but arbitrarily selected) solution to (2) for position values 𝑥(𝑡0 ) near the nominal ones. Thus, if the particular solution is such that the position value 𝑥(𝑡0 ) at time 𝑡0 satisfies 𝑥(𝑡0 ) ∈ 𝑈0 + {𝑥0 } then 𝑣 𝑔 (is well defined at 𝑥(𝑡0 ) and) represents the error, or deviation, from required velocity for the particular solution. From Prop. III.1 it follows that zero error 𝑣 𝑔 = 0 holds along a nominal trajectory and from (47) and (51) it follows more generally that 𝑣 𝑔 at time 𝑡0 has the following properties 𝜕𝑣 𝑔 (𝑥, 𝑡0 ) 𝜕𝑡 𝜕𝑣 𝑔 (𝑥, 𝑡0 ) 𝜕𝑥

= =

𝜕𝑣 𝑟 (𝑥, 𝑡0 ) − 𝑔 − 𝑎(𝑥, 𝑣(𝑡0 )) − 𝑢(𝑡0 ), 𝜕𝑡 𝜕𝑣 𝑟 (𝑥, 𝑡0 ) = 𝑄(𝑥, 𝑡0 ), 𝜕𝑥

𝑥 ∈ 𝑈0 + {𝑥0 }.

By continuity of the solution to (2) with respect to 𝑢 it follows that there is a subset 𝒰0 of the class of continuous controls (containing the control present in the particular solution above) and a time interval ℐ0 containing 𝑡0 such thatr 𝑑 𝜕𝑣 𝑟 (𝑥(𝑡), 𝑡) 𝜕𝑣 𝑟 (𝑥(𝑡), 𝑡) 𝑣 𝑔 (𝑥(𝑡), 𝑡) = + 𝑣(𝑡) − 𝑔 − 𝑎(𝑥(𝑡), 𝑣(𝑡)) − 𝑢(𝑡) = 𝑑𝑡 𝜕𝑡 𝜕𝑥 )︀ 𝜕𝑣 𝑟 (𝑥(𝑡), 𝑡) 𝜕𝑣 𝑟 (𝑥(𝑡), 𝑡) (︀ + 𝑣 𝑟 (𝑥(𝑡), 𝑡) − 𝑣 𝑔 (𝑥(𝑡), 𝑡) − 𝑔 − 𝑎(𝑥(𝑡), 𝑣(𝑡)) − 𝑢(𝑡) = 𝜕𝑡 𝜕𝑥 ˜ (𝑥(𝑡), 𝑡), 𝑡 ∈ ℐ0 , 𝑢 ∈ 𝒰0 , − 𝑄(𝑥(𝑡), 𝑡)𝑣 𝑔 (𝑥(𝑡), 𝑡) − 𝑎 r Note

that we cannot assume that (41) holds here.

15 of 24 American Institute of Aeronautics and Astronautics

(52)

˜ is given by where we have used (44), (47), (51), and 𝑎 (︀ )︀ ˜ (𝑥, 𝑡0 ) = 𝑎(𝑥, 𝑣(𝑡0 )) − 𝑎 𝑥, 𝑣 𝑟 (𝑥, 𝑡0 ) + 𝑢(𝑡0 ), 𝑎

𝑥 ∈ 𝑈0 + {𝑥0 }.

(53)

Downloaded by John Robinson on January 26, 2015 | http://arc.aiaa.org | DOI: 10.2514/6.2015-0087

The relation (52) is the fundamental relation for VTG based guidance. By continuity of the solutions to (2) with respect to both initial values and controls it follows that there exists a neighborhood 𝒩 of the initial condition (𝑥𝐺 , 𝑣 𝐺 ) to the firing solution and a neighborhood 𝒰¯0 of the zero control such that when the solution 𝑡 ↦→ (𝑥(𝑡), 𝑣(𝑡)) starts with initial conditions in 𝒩 and has controls in 𝒰¯0 then the relation (52) holds for ℐ0 = ℐ and 𝒰0 = 𝒰¯0 . ˜ in (53) as Using a simple perturbation expansion we can express 𝑎 (︀ )︀ 𝜕𝑎 𝑥, 𝑣 𝑟 (𝑥, 𝑡0 ) ˜ (𝑥, 𝑡0 ) = − 𝑣 𝑔 (𝑥, 𝑡0 ) + 𝑢(𝑡0 ) + 𝑜(‖𝑣 𝑔 (𝑥, 𝑡0 )‖). (54) 𝑎 𝜕𝑣 For the Jacobian on the right we can apply the drag-only approximation in sec. II.E.1 and use (14). If we do this and use similar arguments about insensitivity to position and velocity as those leading to (49) we arrive at the approximation (︀ )︀ )︀ 𝜕𝑎 𝑥, 𝑣 𝑟 (𝑥, 𝑡0 ) 𝜌(𝑥0 )𝑆𝐶𝐷 ‖𝑣 0 ‖ (︀ 1 𝑣 0 𝑣 𝑇0 , 𝑥 ∈ 𝑈0 + {𝑥0 }. =− 𝐼 3×3 + (55) 2 𝜕𝑣 2𝑚 ‖𝑣 0 ‖ Connection with ZEPD techniques. We continue to consider the situation in sec. III.B.2 (where we assume that 𝑥(𝑡0 ) ∈ 𝑈0 + {𝑥0 }) and recall the definition of the zero effort terminal deviation 𝜁 0 (𝑡𝐹 ) in (26) (where 𝑡𝐹 is an arbitrary terminal time). Using a simple perturbation expansion we can express the ZEPD in (31) at 𝑡𝐹 = 𝑡𝑇 as 𝜁 0,𝑥 (𝑡𝑇 ) = [Φ𝑥𝐺 ((𝑥(𝑡), 𝑣(𝑡)), 𝑡𝑇 − 𝑡]1:3 − 𝑥𝑇 = [Φ𝑥𝐺 ((𝑥(𝑡), 𝑣(𝑡)), 𝑡𝑇 − 𝑡]1:3 − [Φ𝑥𝐺 ((𝑥(𝑡), 𝑣 𝑟 (𝑥(𝑡), 𝑡)), 𝑡𝑇 − 𝑡)]1:3 = −

𝜕[Φ𝑥𝐺 ((𝑥(𝑡), 𝑣 𝑟 (𝑥(𝑡), 𝑡)), 𝑡𝑇 − 𝑡)]1:3 𝑣 𝑔 (𝑥(𝑡), 𝑡) + 𝑜(‖𝑣 𝑔 (𝑥(𝑡), 𝑡))‖), 𝜕𝑣

𝑡 ∈ ℐ0 , (56)

where we have used (45), the definition (51) and ℐ0 is some time interval around 𝑡0 . This is the fundamental relation between ZEPD and VTG. By the assumption in (40), the Jacobian on the right in (56) is nonsingular so if we neglect the higher order term in (56) there is locally a one-to-one relation between 𝜁 0,𝑥 (𝑡𝑇 ) and 𝑣 𝑔 . The fact that the control 𝑢 does not enter on the right hand side of (56) manifests the interpretation of the VTG as the difference between the required velocity for free flight to the the target and the current velocity. If we recall (16) we have for small values of 𝑣 𝑔 and time-to-go the approximation 𝜁 0,𝑥 (𝑡𝑇 ) = −(𝑡𝑇 − 𝑡)𝑣 𝑔 (𝑥(𝑡), 𝑡) +

)︀ (𝑡𝑇 − 𝑡)2 𝜌(𝑥(𝑡))𝑆𝐶𝐷 ‖𝑣(𝑡)‖ (︀ 1 𝐼 3×3 + 𝑣(𝑡)𝑣(𝑡)𝑇 𝑣 𝑔 (𝑥(𝑡), 𝑡). (57) 2 2 2𝑚 ‖𝑣(𝑡)‖

In the vacuum trajectory case 𝑎 = 0 we obtain the simple approximations 𝜁 0,𝑥 (𝑡𝑇 ) = −(𝑡𝑇 − 𝑡)𝑣 𝑔 (𝑥(𝑡), 𝑡) = 𝑄(𝑥(𝑡), 𝑡)−1 𝑣 𝑔 (𝑥(𝑡), 𝑡).

(58)

This clearly illustrates the relation between the two basic forms of error that can be used to drive the guidance loop. In particular, it shows that the VTG is not a “derivative” of the ZEPD but is in fact an error quantity on the same derivative “level” in the sense that if one of these quantities is zero at some time instant, the other is too. Therefore, there is no difference in formulating control problems in terms of ZEPD or VTG, but the latter should offer some advantage due to the simple and intuitive dynamic relation (52). Control problem for the velocity to be gained. Due to the equivalence between the ZEPD and the VTG expressed by (56)–(58) the guidance problem can be formulated around either of these two variables. The simple form of the dynamics relation (52) for the VTG 𝑣 𝑔 makes it is reasonable to take this as a starting point, representing the open loop dynamics for a state feedback control problem. s When

the time-to-go becomes small (or the ballistic coefficient is large) there is little difference between (57) and (58).

16 of 24 American Institute of Aeronautics and Astronautics

From (52) and (54) obtain an approximation of the open loop dynamics for 𝑣 𝑔 as (︀ )︀ (︁ 𝜕𝑎 𝑥(𝑡), 𝑣 𝑟 (𝑥(𝑡), 𝑡) )︁ 𝑑 𝑣 𝑔 (𝑥(𝑡), 𝑡) − 𝑢(𝑡). 𝑣 𝑔 (𝑥(𝑡), 𝑡) = − 𝑄(𝑥(𝑡), 𝑡) − 𝑑𝑡 𝜕𝑣

(59)

The model in (59) is valid near the nominal trajectory and for small controls, in the same sense as (52). If we make a change of coordinates on the control as (︀ )︀ (︁ 𝜕𝑎 𝑥(𝑡), 𝑣 𝑟 (𝑥(𝑡), 𝑡) )︁ 𝜈(𝑡) = 𝑢(𝑡) + 𝑄(𝑥(𝑡), 𝑡) − 𝑄0 (𝑥(𝑡), 𝑡) − 𝑣 𝑔 (𝑥(𝑡), 𝑡) (60) 𝜕𝑣

Downloaded by John Robinson on January 26, 2015 | http://arc.aiaa.org | DOI: 10.2514/6.2015-0087

we can write (59) as 𝑑 𝑣 𝑔 (𝑥(𝑡), 𝑡) = −𝑄0 (𝑥(𝑡), 𝑡)𝑣 𝑔 (𝑥(𝑡), 𝑡) − 𝜈(𝑡), (61) 𝑑𝑡 where 𝑄0 is the 𝑄-matrix in (50). This is the standard form for the VTG dynamics in vacuum (and constant gravity) which is a (completely) controllable system (for 𝑡 < 𝑡𝑇 ). An approximation of the matrix factor on the right in (60) can be obtained from (49) and (55) ast (︀ )︀ )︀ 𝜕𝑎 𝑥(𝑡), 𝑣 𝑟 (𝑥(𝑡), 𝑡) 1 𝜌(𝑥(𝑡))𝑆𝐶𝐷 ‖𝑣(𝑡)‖ (︀ 𝐼 3×3 + (62) = 𝑣(𝑡)𝑣(𝑡)𝑇 . 𝑄(𝑥(𝑡), 𝑡) − 𝑄0 (𝑥(𝑡), 𝑡) − 𝜕𝑣 4𝑚 ‖𝑣(𝑡)‖2 In analogy with a remark made in connection with the CPA in (32) we note that a nonzero component of VTG in the direction to the target is essentially is equivalent to a change in time of arrival, when the time-to-go is small. Therefore, it is motivated to make a change of coordinates in (61) to a frame where the 𝑥-axis is parallel with the predicted terminal velocity vector 𝑣 𝑇 , and the two other axes provide coordinates in [𝑣 𝑇 ]⊥ . The VTG is then projected onto the plane [𝑣 𝑇 ]⊥ . This is asymptotically, for small values of the time-to-go,u equivalent (cf. (35) and (57)) to substituting 𝜁 0,𝑥 (𝑡𝐶𝑃 𝐴 ) for 𝜁 0,𝑥 (𝑡𝑇 ) in (56). With such a procedure, the VTG dynamics in (61) effectively decouples into two scalar equations of the form 𝑥(𝑡) ˙ = 𝑄0 (𝑡)𝑥(𝑡) + 𝑢(𝑡), 𝑄0 =

1 , 𝑡𝑇 − 𝑡

𝑡 ∈ [𝑡0 , 𝑡𝑇 ),

(63)

for the lateral and longitudinal directions in [𝑣 𝑇 ]⊥ (since the axial VTG component is ignored). Linear quadratic optimal control. In the guided projectile application a significant part of the deviations from a nominal trajectory normally stem from errors in initial velocity and aiming of the gun. Moreover, early control actions have large “leverage”. Therefore, it is natural to seek a control formulation that time-weights the error in such a way that a significant portion of the control action is performed early. A straightforward way to accomplish this is to apply (finite horizon) linear quadratic (LQ) control. A standard LQ control problem32 for the system in (63) can be formulated as to seek a control 𝑢( · ) ∈ 𝒰, where 𝒰 is e.g. the class of piecewise continuous controls on [𝑡0 , 𝑡1 ], which minimizes the performance functional 𝒱(𝑥(𝑡0 ),𝑡0 ) : 𝒰 → [0, ∞) given byv 2

∫︁

𝑡1

𝒱(𝑥(𝑡0 ),𝑡0 ) (𝑢( · )) = 𝑞1 𝑥 (𝑡1 ) +

(︀

𝑡0

)︀ 𝑞𝑥 𝑥2 (𝑡) + 𝑢2 (𝑡) 𝑑𝑡, 2 (𝑡𝑇 − 𝑡)

(64)

where 𝑡1 < 𝑡𝑇 is a time shortly before the nominal terminal time 𝑡𝑇 . The constants 𝑞1 , 𝑞𝑥 > 0 are penalty weights to be specified below. The form on the penalty term for the state in the integrand in (64) is chosen so that errors near the right end point of [𝑡0 , 𝑡1 ] incur a higher cost. The penalty on control effort is uniform, since the actuators work under similar conditions during all of the guided phase. A minimizing (continuous) control 𝑢 ˆ for the problem in (64) exists in this case and is of the form 𝑢 ˆ(𝑡) = −𝑝(𝑡)𝑥(𝑡),

(65)

to the size estimates in Ref. 25 the quantity in (62) is of the order 10−2 for an 155mm artillery scenario. larger values of the time-to-go it can be meaningful to instead use the current velocity 𝑣(𝑡) and the plane [𝑣(𝑡)]⊥ to construct a (rotating) frame in which the equation (61) can be expressed. v For a related LQ problem see p. 154 in Ref. 33. t According

u For

17 of 24 American Institute of Aeronautics and Astronautics

where 𝑝 is the (unique) positive solution to the Riccati equation − 𝑝(𝑡) ˙ =

2 𝑞𝑥 , 𝑝(𝑡) − 𝑝2 (𝑡) + 𝑡𝑇 − 𝑡 (𝑡𝑇 − 𝑡)2

𝑝(𝑡1 ) = 𝑞1 ,

𝑡 ∈ [𝑡0 , 𝑡1 ],

(66)

It is easy to verify that if 𝑞1 and 𝑞𝑥 satisfy √︀ 3(1 + 1 + 4𝑞𝑥 /9) 𝑞1 = 2(𝑡𝑇 − 𝑡1 )

(67)

then the (unique, positive) solution to the Riccati equation (66) is given by

Downloaded by John Robinson on January 26, 2015 | http://arc.aiaa.org | DOI: 10.2514/6.2015-0087

𝑝(𝑡) = 𝑞1

(𝑡𝑇 − 𝑡1 ) . (𝑡𝑇 − 𝑡)

(68)

It is natural to consider 𝑞1 as the main parameter to vary when “tuning” the controller (so that 𝑞𝑥 is determined by (67)) since it represents penalty on a variable (a component of the velocity to be gained near the nominal terminal time) which has a clear physical meaning (cf. e.g. (57)). With the optimal control 𝑢 ˆ given by (65), (67) and (68) the closed loop guidance dynamics will bew 𝑥(𝑡) ˙ = (𝑄0 (𝑡) − 𝑝(𝑡))𝑥(𝑡) =

1 − 𝑞1 (𝑡𝑇 − 𝑡1 ) 𝑥(𝑡). (𝑡𝑇 − 𝑡)

This is a simple linear time varying first order system (which, for each fixed 𝑡, represents stable dynamics due to (67)). III.B.3.

Perturbation based velocity to be gained guidance algorithm

An algorithm for velocity to be gained guidance based on the perturbation expressions for zero effort position deviation can now be given as follows. Before flight, compute a nominal trajectory (23) using any kind of high fidelity method. Then, at each guidance update instant during flight, do: (1) Perform the first two steps of the algorithm for predictive guidance in section III.A.3 to obtain an approximation to the CPA 𝜁 0,𝑥 (𝑡𝐶𝑃 𝐴 ) in (32). (2) Compute the velocity to be gained in the plane [𝑣 𝑇 ]⊥ perpendicular the nominal terminal velocity 𝑣 𝑇 . This is done by substituting 𝜁 0,𝑥 (𝑡𝐶𝑃 𝐴 ) for 𝜁 0,𝑥 (𝑡𝑇 ) in (57), and solving for 𝑣 𝑔 . Extract the lateral and longitudinal components of 𝑣 𝑔 (based on coordinates in the plane [𝑣 𝑇 ]⊥ ). (3) Compute lateral and longitudinal acceleration components in the plane [𝑣 𝑇 ]⊥ for the control 𝜈 in (61) using the corresponding components of 𝑣 𝑔 and the generic formula (65) (where 𝑥 plays the role of a component of 𝑣 𝑔 and −𝑢 is the corresponding acceleration command). Map the control values for 𝜈 back to the physical control values 𝑢 using (60). This formulation yields a closed loop guidance for continuous time (i.e. short sampling interval) operation.

IV.

Simulations

We shall now illustrate the results above with some simulations. The projectile dynamics are described with a 6DOF model and the differential equations are integrated using a fourth-order Runge-Kutta ODE solver with a time step of 0.01s. (The same equations of motion and ODE solver are used for computation of a nominal unguided trajectory and for all the unguided and guided perturbed trajectories.) The Earth is considered as non rotating. w This assumes that the instantaneous guidance loop time constant −1/(𝑄 (𝑡) − 𝑝(𝑡)) never becomes too small. If 𝑞 , 𝑡 are 0 1 1 chosen such that −1/(𝑄0 (𝑡) − 𝑝(𝑡)) can approach the effective time constant in the projectile dynamics then a static autopilot together with a simple command filter may not be enough and it may be necessary to have a dynamics synthesizing autopilot.

18 of 24 American Institute of Aeronautics and Astronautics

IV.A.

Shell and firing scenario

The firing scenario is the same as in Ref. 25 with gun elevation angle 30∘ , initial nominal (muzzle) projectile velocity 684.3m/s and initial spin rate 1102.6rad/s. The nominal downrange distance traveled is approximately 16km and the (rightwards) crossrange drift 296m. A projectile model based on a standard M107 artillery shell is used with mass properties, dimensions and aerodynamic coefficients taken from Ref. 34. Data for the shell is summarized in table 1. Table 1: Shell properties.

Downloaded by John Robinson on January 26, 2015 | http://arc.aiaa.org | DOI: 10.2514/6.2015-0087

Description Total mass Caliber Axial moment of inertia Lateral moment of inertia

IV.B.

Value 43.0kg 0.155m 0.144kg m2 1.216kg m2

Fuze and guidance properties

The fuze is modeled as a generic 2D course correction fuze which is roll stabilized with respect to the horizontal plane. It can generate forces in the two directions perpendicular to the main axis, representing the effect of normal forces from deflectable canards mounted near the nose. The moment lever arm with respect to the center of mass is set to 0.4m. A simple linear second order model for the actuator dynamics is used to represent the response of the actuators to the commands from the autopilot (see below). The data for the actuators is given in table 2. The data for the available control force from canards is taken from Ref. 6. A limit of 60N for the canard control force for both the horizontal and vertical directions is used in all simulations except those relating to the influence of the control force limit. The value 60N is within the range of values obtained for the configurations considered in Ref. 6 and is well within the achievable range with respect to packaging and stability requirements.6 The guidance commands are computed using the steps outlined in section III.B.3 above and the guidance loop operates in continuous time with parameters 𝑞1 = 8 and 𝑡1 = 0.5s (cf. (64)–(68)). Guidance is initiated after 15s of flight and operates until impact but with 𝑝(𝑡) in (65) being “frozen” during the last 0.5s. A simple autopilot, described in Ref. 25, 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 uses a perfect model of the projectile and ideal measurements, and it does not include any dynamics synthesis. Instead, a simple second order low-pass filter, described in Ref. 25 and summarized in table 2, is inserted to act on the guidance commands before these are passed to the autopilot. The guidance filter parameters are the same as the one used in Ref. 25 except that the guidance filter natural frequency and commanded acceleration limit are both lower here, see table 2. Table 2: Guidance filter parameters (left) and actuator parameters (right). Description Commanded acc. limit Commanded acc. rate limit Guidance filter natural freq. (𝜔𝑔 ) Guidance filter damping (𝜁𝑔 )

IV.C.

Value ±0.5 · 9.81m/s2 ±0.02 · 9.81m/s2 /s 0.5 · 2𝜋rad/s 0.8

Description Canard control force limit Canard ctrl. force moment arm Actuator natural frequency (𝜔𝑎 ) Actuator damping (𝜁𝑎 )

Value ±60N 0.4m 4 · 2𝜋rad/s 0.8

Dispersion analysis

Let 𝑥𝑘 (𝑡𝑖 ) denote the coordinates in the fire control frame 𝐹 of the 𝑘th run at the time of impact 𝑡𝑖 with the ground plane. The miss distance 𝑑𝑘 for the 𝑘th simulation run is then defined as 𝑑𝑘 = ‖𝑥𝑘 (𝑡𝑖 ) − 𝑥𝑇 ‖, where 𝑥𝑇 denotes the coordinates in 𝐹 of the nominal ground impact point. (For the downrange and crossrange coordinates, mean and standard deviation are defined in the standard way with error with sign.)

19 of 24 American Institute of Aeronautics and Astronautics

The dispersion is analyzed using Monte Carlo simulations with Gaussian zero mean uncorrelated errors introduced to the initial velocity, elevation and azimuth. The standard deviation for the initial velocity error is set to 3m/s and the standard deviations for the angle errors are both set to 0.1∘ . In table 3 the dispersion statistics for an unguided shell is presented for different combinations of perturbations in initial velocity, elevation and azimuth angles. The corresponding statistics for a guided shell is presented in table 4. The guided shell is launched at an aimpoint approximately 40m beyond the desired impact point 𝑥𝑇 (cf. section IV.F) but the nominal trajectory used for guidance is one with impact point at the exact target location 𝑥𝑇 .

Downloaded by John Robinson on January 26, 2015 | http://arc.aiaa.org | DOI: 10.2514/6.2015-0087

Table 3: Statistics from 1000 Monte Carlo simulations of the unguided projectile with perturbations on the initial elevation, azimuth and velocity. Perturbation in Elevation Mean Std. Dev. 17.686 13.401 0.406 22.139 0.034 1.543

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

one parameter Velocity Mean Std. Dev. 59.943 45.439 -3.619 75.119 -0.108 2.347

Perturbation in several parameters Elev. and vel. El., vel. and az. Mean Std. Dev. Mean Std. Dev. 63.113 46.848 70.820 43.979 -3.213 78.509 -3.256 78.517 -0.074 2.821 0.956 27.897

Table 4: Statistics from 1000 Monte Carlo simulations of the projectile with course correcting fuze using velocity to be gained guidance and perturbations added on the initial elevation, azimuth and velocity.

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

Perturbation in Elevation Mean Std. Dev. 0.492 0.357 -0.483 0.364 -0.057 0.027

one parameter Velocity Mean Std. Dev. 1.065 4.291 -0.268 4.406 -0.040 0.248

Perturbation in several parameters Elev. and vel. El., vel. and az. Mean Std. Dev. Mean Std. Dev. 1.234 4.836 1.517 5.652 -0.341 4.971 -0.242 5.835 -0.040 0.281 -0.030 0.380

In figure 4a and 4b the dispersion is presented in form of histograms. These two figures represent the dispersion for perturbations corresponding to the case presented in the rightmost columns of table 3 and 4 (including perturbations in elevation, velocity and azimuth). Notice that while the dispersion for the guided shell is very low, there are some outliers. These can be traced back to control force saturation (see figure 8c). In figure 7 these outliers are more clearly visible.

(a) Unguided shell.

(b) Guided shell.

Figure 4: Error distribution of the impact point based on 1000 simulated trajectories.

20 of 24 American Institute of Aeronautics and Astronautics

Downloaded by John Robinson on January 26, 2015 | http://arc.aiaa.org | DOI: 10.2514/6.2015-0087

(a) Unguided shell.

(b) Guided shell.

Figure 5: Influence of separate initial errors on an unguided and a guided shell as a function of elevation.

IV.D.

Dispersion over elevation

In order to illustrate the performance of the proposed guidance strategy for varying firing scenario parameters, Monte Carlo simulations were conducted for different nominal elevations. The influence of errors in initial velocity, elevation and azimuth are analyzed separately. The errors in launch parameters are also here Gaussian, zero mean and uncorrelated with standard deviation for the velocity and angle errors of 3m/s and 0.1∘ , respectively. As can be seen in figure 5b, the guidance and control strategy is capable of significantly reducing the dispersion as long as sufficient flight time is achieved. For an elevation of 10∘ the nominal flight time is just above 20s and since the guidance does not commence until after 15𝑠 very small improvements are obtained compared to the unguided case. The aimpoint is here the same as the nominal impact point 𝑥𝑇 for the different elevations. The results in figure 5b were generated with a 𝑞1 = 6 which produced more stable results over all elevations but slightly less accurate for the baseline of 30∘ elevation. If required, the guidance parameter 𝑞1 can be scheduled prior to launch. IV.E.

Limitations in control force

As stated above, a control force limit of 60N was assessed to be a realistic value for the selected configuration and used as a baseline for the presented results. In figure 6 the control force is limited to evaluate the impact

Figure 6: Effect of control force limit on downrange standard deviation.

21 of 24 American Institute of Aeronautics and Astronautics

Downloaded by John Robinson on January 26, 2015 | http://arc.aiaa.org | DOI: 10.2514/6.2015-0087

Figure 7: Result of aimpoint on dispersion with resulting 2𝜎 error covariances. The bottom panel represents aiming at the nominal impact point and the top panel aiming approximately 40m beyond.

on downrange errors. The same scenario as described in section IV.A is used here. It can be seen that even lower limits still give an significant improvement in downrange dispersion compared to an unguided shell, and for control forces above 80N no apparent improvement is seen (for this scenario with a 30∘ nominal elevation). IV.F.

Aimpoint

The guided shell has an inherent difficulty to gain range due to the limited lift force possible to generate for the projectile. This results in some outliers in impact point when aiming at the nominal impact point, as illustrated in bottom of figure 7, mainly for the cases when the initial velocity is significantly low (around 8 − 10m/s below nominal). This can be managed by aiming slightly beyond the desired impact point, which is shown in the top of figure 7 (here approximately 40m beyond, corresponding to an increase in elevation by 0.2∘ ). IV.G.

Guidance loop behavior

In figure 8 the control force time history and aerodynamic angles are presented for the baseline simulation case presented in section IV.A. From the 1000 Monte Carlo runs every 20th was chosen to avoid cluttering. The total angle of attack are kept within an outer bound of just below 4∘ , indicating a stable flight. Notice that for a few of the trajectories the control force reaches the 60N limit.

V.

Conclusion

The proposed velocity to be gained guidance has been shown to perform very well for errors of moderate size, corresponding to the type of errors that can occur in real firing situations. However, the outlier cases with extreme velocity errors present difficulties for this combination of guidance and fuze concept. One of the reasons for this is the limited control authority available and the limited additional lift that can be provided by the shell in case the initial velocity is very different compared to the nominal. Another reason is that the outlier cases generate large control errors and corresponding large commanded accelerations, which may excite the projectile dynamics and make the implicit assumption of time scale separation between the guidance and control (and plant) invalid. For these situations a more sophisticated autopilot could improve the performance. As for the implementation of the algorithm it can be noted that while the derivation of the guidance equations is somewhat lengthy, the result is very easy to implement, is of low complexity, and does not require sophisticated sensors on the fuze. In fact, since an advanced autopilot is not required, acceleration measurements need only be updated on the (relatively slow) time scale that the guidance operates on.

22 of 24 American Institute of Aeronautics and Astronautics

Downloaded by John Robinson on January 26, 2015 | http://arc.aiaa.org | DOI: 10.2514/6.2015-0087

(a) Velocity to be gained, 𝑁 frame.

(b) Total angle of attack.

(c) Control force, 𝑁 frame.

(d) Aerodynamic angles.

Figure 8: Time history of guidance and control variables. In (a) the velocity to be gained variable is illustrated for one arbitrarily selected trajectory (the 900th) and in (b), (c) and (d) it is composed of every 20th trajectory from the original set of 1000.

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

References 1 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. 2 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. 3 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. 4 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. 5 Fresconi, F. and Plostins, P., “Control Mechanism Strategies for Spin-Stabilized Projectiles,” Tech. Rep. ARL-TR-4611, US Army Research Laboratory, September 2008. 6 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). 7 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. 8 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.

23 of 24 American Institute of Aeronautics and Astronautics

Downloaded by John Robinson on January 26, 2015 | http://arc.aiaa.org | DOI: 10.2514/6.2015-0087

9 Ollerenshaw, D. and Costello, M., “Model Predictive Control of a Direct Fire Projectile Equipped with Canards,” Proc. AIAA Guidance, Navigation, and Control Conference, San Francisco, CA, August 2005, AIAA paper 2005-5818. 10 Gagnon, E. and Lauzon, M., “Maneuverability Analysis of the Conventional 155 mm Gunnery Projectile,” Proc. AIAA Guidance, Navigation and Control Conference, Hilton Head, SC, August 2007, AIAA paper 2007-6784. 11 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, Honolulu, HI, August 2008, AIAA paper 2008-6997. 12 Calise, A. and El-Shirbiny, H., “An Analysis of Aerodynamic Control for Direct Fire Spinning Projectiles,” Proc. AIAA Guidance, Navigation and Control Conference, Montreal, Canada, August 2001, AIAA paper 2001-4217. 13 Gross, M. and Costello, M., “Impact Point Model Predictive Control of a Spin-Stabilized Projectile with Instability Protection,” Proc. AIAA Atmospheric Flight Mechanics Conference, Boston, MA, August 2013, AIAA paper 2013-4509. 14 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, Portland, OR, August 2011, AIAA paper 2011-6247. 15 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. 16 S` eve, F., Theodoulis, S., Wernert, P., Zasadzinski, M., and Boutayeb, M., “Pitch/Yaw Channels Control Design for a 155mm Projectile with Rotating Canards, using a 𝐻∞ Loop-Shaping Design Procedure,” Proc. AIAA Guidance, Navigation and Control Conference, National Harbor, MD, January 2014, AIAA paper 2014-1474. 17 Salman, M. and Chang, B., “Active Coning Compensation for Control of Spinning Flying Vehicles,” Proc. IEEE Intl. Conference on Control Applications, Yokohama, Japan, September 2010. 18 Laning, J. and Battin, R., “Theoretical Principle for a Class of Inertial Guidance Computers for Ballistic Missiles,” Tech. Rep. R-125, MIT Instrumentation Laboratory, Cambridge, MA, June 1956. 19 Battin, R., “Some New Analytical Results of Use in the Calculation of Ballistic Missile Trajectories,” Tech. Rep. R-181, MIT Instrumentation Laboratory, Cambridge, MA, July 1958. 20 Townsend, G., Abbott, A., and Palmer, R., “Guidance, flight mechanics and trajectory optimization,” Tech. Rep. NASACR-1007, NASA, April 1968. 21 Battin, R., An Introduction to the Mathematics and Methods of Astrodynamics, AIAA, Reston, VA, revised ed., 1999. 22 Rinalducci, A., “A Launch Strategy fo High-Inclination Orbit Missions from the San Marco Equatorial Range,” Acta Astronautica, Vol. 41, 1997, pp. 559–568. 23 Bhat, M. and Shrivastava, S., “An Optimal Q-Guidance Scheme for Satellite Launch Vehicles,” J. Guidance, Vol. 10, No. 1, Jan.-Feb. 1987, pp. 53–60. 24 Burns, S. and Scherock, J., “Lambert Guidance Routine Designed to Match Position and Velicity of Ballistic Target,” J. Guidance, Control, and Dynamics, Vol. 27, No. 6, November-December 2004, pp. 989–996. 25 Robinson, J. W. and Str¨ omb¨ ack, P., “Perturbation Based Guidance for a Generic 2D Course Correcting Fuze,” Proc. AIAA Guidance, Navigation and Control Conference, Boston, MA, August 2013, AIAA paper 2013-5114. 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 Teschl, G., Ordinary Differential Equations and Dynamical Systems, Vol. 140 of Graduate Studies in Mathematics, American Mathematical Society, providence, RI, 2012. 28 Khalil, H., Nonlinear Systems, Prentice-Hall, Upper-Saddle River, NJ, 3rd ed., 2002. 29 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). 30 Brauer, F., “Perturbations of Nonlinear Systems of Differential Equations. II,” J. Math. Anal. Appl., Vol. 17, 1967, pp. 418–434. 31 Sarture, C., “Guidance Schemes for Rocket Vehicles,” Tech. Rep. STL/TN-60-0000-GR186, TRW Space Technology Laboratories, Los Angeles, CA, October 1960. 32 Anderson, B. and Moore, J., Linear Optimal Control, Prentice-Hall, Englewood Cillfs, NJ, 1971. 33 Bryson, A. and Ho, Y.-C., Applied Optimal Control, Taylor-Francis, New York, NY, 1975, Revised printing. 34 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.

24 of 24 American Institute of Aeronautics and Astronautics