Hitting a Moving Target: The Missile Guidance System

• Download Guided_Missile_Simulator source files and demo - 130.9 KB
• Download Time_Graph_GUI source files and demo - 99.5 KB
• Download Simple_Guidance_System_UI source files and demo - 106.4 KB

Overview of Downloads

The simulator itself is a console window. The kernel beep Doppler effect output and the HTML copy of the console window saved to the desktop can be turned on or off by adjusting the values of the Boolean constants at the top of Program.cs.

The GUI helps visualize the way adjusting variables alters the time it will take to intercept the target. Remember, where the red/blue line crosses the x axis is an intercept time and when the green/purple line crosses the x axis is a valid closest approach, which if less than the radius of tolerance, is a glancing blow, and not a miss.

The simple UI provides a sample of a simple application which uses my guidance system in a graphical interface. The concept is quite simple: you are the little red spaceship that starts out in the middle of the screen, and the green spaceships are trying to hit you with their white missiles. Note that anything going too far in any direction will wrap around to the other side to prevent it from flying off screen. Also note that your ship and the green ships have a limit on their velocities to prevent the situation from getting out of hand. Additionally, your ship will turn orange when you thrust forward or backward, but not when that thrust is thrown out because it would exceed the velocity limit.

Introduction

O.K., here I go with my first article. I hope you like it; I've put my heart and soul into it. Comments are more than welcome, and you can count on me to read them and respond to any and all questions, comments, or concerns.

I became fascinated with guided missiles back in August of 2006 when I attended a game programming workshop. I came up with a crude way of leading a target because I did not like that the enemies in my game shot at where the player was, rather than where the player would be when the bullet got to him. My solution worked somewhat, but it was strictly limited to two dimensions and it produced a rough estimate of the correct firing angle through the use of a loop; I knew this was algorithmically inefficient, but I had a 72 hour time-line, and what I threw together would have to suffice.

After that time I researched the problem on my own. Time was limited because I was starting my senior year in high school, but I was fortunate enough to be able to use this topic as my senior project, a project with a self-assigned topic that every student must, it is hoped, learn something through. I did not know how to solve this complex mathematical problem, and I was able to get some time set aside for my research as a result.

This project summarizes my findings. I have developed a system for intercepting a target with a position, velocity, and acceleration with a missile (fired from a shooter with a given position and velocity-- shooter acceleration does not effect the missile) which has a given acceleration for a given duration of time (think of it as fuel, which anticipates no acceleration after being exhausted), a given initial or "muzzle" velocity, and a given initial radial distance (as in the length of the launcher which the missile appears at the end of when fired, with said initial velocity). The results of this calculation come in the form of a Boolean return indicating the ability for the specified missile to hit the specified target from the specified shooter, a list of unit vectors (a 2D double array) in which the first index is the intercept direction (there may be up to four possibilities!) listed in ascending order (fastest intercept is first in the array) and the second index (0, 1, or 2) indicates X, Y, or Z respectively, and, optionally, how long it will take each of these vectors to impact the target. My findings are free for public use, and incredibly useful for game programming, especially for developing good AI for a space shooter, among other things. (Please remember to give Christopher Kampas credit (please, call me Jay)) if you do choose to use my guidance system DLL or any other component of my project for your own purposes, in fact I highly encourage it, as the purpose of this project was to make life easier for a game developer, most certainly not just me.)

Background

Relax, while I have used basic calculus gospel and the quartic formula blended together in complex algebra, you, the developer, are not expected to know the first thing about these topics, though it certainly wouldn't hurt if you do (conceptually, understanding polar graphing couldn't hurt either). I will provide here the basic mathematical foundation required to understand what my project does; when we're done here, I hope, the mathematics will be completely demystified. Essentially I solve for the amount of time it will take to hit the target, and then simply figure out where it will be after that amount of time has passed and aim in that direction.

First of all, if I want to determine the path of interception relative to the shooter (with acceleration disregarded because it does not effect the initial state of the missile as velocity and position do), by subtracting the shooter's position and velocity from the target's respective values the result is a scenario that will have the same intercept vector(s) as our original scenario, but firing from the origin at these new target values instead. We have already completely eliminated the shooter's parameters from the problem.

Next, handed down from calculus, are two equations that help us describe the position of an object as a function of time. Firstly, the velocity of an object after given time T has passed is equal to the acceleration times T, plus the initial velocity (Vt = At + V0). Second, the position of an object after given time T has passed is equal to one half of the acceleration times T squared, plus the initial velocity times T, plus the initial velocity (Pt = (1/2)At2 + V0t + P0). Specifically why these equations work is somewhat complicated to explain (especially the second one, with it's mysterious 1/2 appearing from seemingly nowhere... but the way a teacher of mine explained it in detail to me, if you think area of a triangle you will be going in the right direction)

Lastly, and defiantly the most obscure piece is the use of the quartic formula. Many of you may know the quadratic formula, which states that where 0 = AX2 + BX + C it follows that X = (-B ± √(B2 - 4AC)) / 2A, a conclusion reached by "completing the square." Well, the quartic formula is far more complicated, and uses the cubic formula extensively, which results in four solutions. If you want to know more, you can Google "quartic formula" and/or look at the "rootsOf" function in my guidance system DLL, but it is far too large and complicated to reasonably provide here. However, what you need to know is that it solves for X where 0 = AX4 + BX3 + CX2 + DX + E.

Using the Code

First of all, to prevent my code from becoming needlessly long and difficult to sift through, my variable names are quite short but still very intuitive. I use three letters to identify most variables in my guidance system. The first letter indicates whom we are speaking of (T for target, S for shooter, and P for projectile). The second letter indicates the vector we are interested in (A for acceleration, V for velocity, and P for position). Finally, the third letter indicates the particular information we want to know about (X, Y, or Z for the dimensional components, or M for magnitude-- perhaps better understood to the unfamiliar as the length, which, by the way, is always equal to √(X2 + Y2 + Z2). It follows that "tay" is understood as "targetAcclerationY" but makes the code far more readable, as calls to "tay" are VERY frequent. Another example is "pvm," or "projectileVelocityMagnitude." There is one special case, "rbt" (remaining burn time, i.e. fuel remaining), also "rot" (radius of tolerance, i.e. the radius of the target plus the radius of the missile) makes appearances outside of the guidance system in my project, and is necessary for determining impacts which are glancing blows, and is also necessary when the missile is out of fuel, as the call to the intercept method will provide 0 for all projectile vector magnitudes, which results in an otherwise erroneous response of missing the target by a micrometer because there is no way to counteract the minor rounding errors occurring in the background of double precision floating point number arithmetic if the target is an infinitesimally small point, i.e. a radius of tolerance 0.

The guidance system of my project is in a DLL which contains the static class Targeting, which contains CalculateAcceleration and CalculateVelocity methods (simply solving a system of equations that given three successive position values, calculate the current acceleration and velocity of that object respectfully, used once per dimension. Targeting also contains two overloads of Intercept, so that you do not have to provide an out array if you do not care how long each vector given in the 2D vector array will take to impact (the order in the arrays corresponds exactly). the one and only private method is the rootsOf method, which solves a quartic for the independent variable (time in our case, X listed in the background section)

An Abstract Look at the Algorithm

For a more in-depth analysis, please read the actual lower-level mathematics is described in the algebraic solution section below. This section aims at conveying the concept abstractly in a simple 2D scenario, and does not discuss the specifics of how the actual results are calculated in my DLL.

Firstly, subtract the positions and velocities of the shooter from the shooter and the target

Position subtraction:

Velocity subtraction:

The next step is to find all occurrences where the missile can be at the same place as the target at the same time (now fired from the origin with no shooter-induced drift at the target's relative position and velocity). Essentially, these are the precise times where the target's location along its parabolic trajectory at some given time intersects the sphere centered at the origin with a radius of the distance that the projectile will have traveled in that particular amount of time (the red target dot and the black circle in the next animation, respectively); there can be as many as four such times, as few as zero, or anywhere in between (the reason why is investigated in the algebraic section), but in practice, due to the nature of the problem, there will simply be zero or two such times in the vast majority of scenarios.

Finally the launch unit vector(s) are calculated (and sorted in order of shortest to longest time in flight). This is done by taking all times calculated in the previous step and determining where the target will be at each of these times. By doing this the missile's guidance system knows where the target will be at every moment in time that it can hit it by traveling in a straight line, so, to calculate the launch unit vectors, each of these coordinates is treated as a vector and normalized. Now all that the caller will have to do is use these vectors to appropriately configure the actual missile object's state and introduce it into their runtime environment (physics engine etc.); the specifics of this environmental configuration can be found in the penultimate paragraph of the algebraic section.

Note that the second stage (the post-rbt "out of fuel" stage is calculated in a very similar way, picking up exactly where the first stage leaves off using a slightly different equation for projectile distance from the origin as a function of time.

As for those so-called "closest" approaches, those equations simply do some calculations to determine the exact points in time when adjusting the value of time used to determine the target's future position to fire at will result in either missing the target by more or missing the target by less, regardless of weather that adjustment was an increase or a decrease (i.e. that is to say that it is the best of times or it is the worst of times, so to speak), and determines if any of these apexes were close enough. Used with the intercept calculations for both stages, these four interlocking equations account for all possible situations and are completely failsafe when used together.

Note that this is not only useful for determining glancing blows off of targets that almost escaped but were too large for their own good, but also for eliminating the adverse side-effects of rounding error when the projectile does not have any acceleration of its own (i.e. pam = 0, which can occur if the missile is nothing more than a constant velocity projectile in and of itself, or, and by far more commonly, because the missile's remaining burn time has reached 0 and it is requesting a trajectory update). To understand this more, experiment with simulation states with these kinds of settings in my GUI download, and take note that the red/blue intercept line will arc towards the x axis and seem to just barely touch it--theoretically it touches at a single, infinitesimally small point; however, due to the rounding error that occurs when doing arithmetic with double-precision floating point numbers will cause the line to either just barely cross it (this is not inherently bad) or it will just barely miss the x axis and begin to pull away from it again without ever touching it (this is exceedingly undesirable, as the guidance system would report that it is unable to intercept, when in fact it will intercept, so long as the target is not infinitesimally small.

Walk-through of the Algebraic Solution

Please note that the above naming conventions will be used here.

Imagine this: After eliminating the shooter's parameters (as in paragraph two of the background section) and applying the equation that determines the position after a given amount of time has passed (Pt = (1/2) * Pam * t2 + Pvm * t + Ppm), it can be said for certain that, so long as I have enough remaining burn time, the distance the missile will be from the origin will be given by this equation; the direction the missile will be fired in is unknown at this time (and will be unknown until the very last step), so while I know how far I will have gone, it could be in any direction at this point. What this means (and why polar graphing will help visualize this, expect diagrams here shortly) is that we can draw a circle in 2D, or a sphere in 3D, etc. which describes where the missile can be after some given amount of time has passed, with the actual location simply being dependant on what direction it was fired. Here is the bottom line: we have an equation that describes the missile's distance from the origin as a function of time.

Next we can use this same formula to describe the target's distance from the origin in every dimension (Xt = (1/2) * Tax * t2 + Tvx * t + Tpx, for instance. the equations for Y and Z are the same with their respective T_y and T_z variables). We can use this with the equation for magnitude to find that Tt = √(Xt2 + Yt2 + Zt2). Note that this is similar to the Pythagorean theorem being done twice as in M = √(√(X2 + Y2)2 + Z2), where the square of a square root simply cancels out.

It follows that when a value for t satisfies the condition Tt = Pt we have found a value for t that can put the missile on the target at the same point at the same time (assuming I fire in the right direction) such that the missile will fly in a straight line during the entire duration of its flight if no evasive action is taken by the target and there is no change in acceleration of either the target or the missile by some outside force; if such a situation would occur, the missile would realign itself as it does every time-step, but this straight line might be slightly bent or the intercept might be made impossible (especially by effective evasive action which uses raw power to escape, i.e. turning away from the missile and opening the throttle up to outrun the missile, assuming the target has enough raw power to achieve this). Ultimately, however, if the acceleration vectors are not changed in any way an intercept is possible, and it will occur after time t has passed from the launch of the missile (assuming we have enough remaining burn time, but I will detail those equations later). Here comes the tough part, solving for t.

While the highest power of t I see when I write out our T and P functions of time is two (i.e. I see time squared, but not anything like time cubed for example), but there is, however, a major problem if we want to try to use the quadratic formula-- the square root in the T equation that results when an intercept is attempted in more than 1D space (in 1D the square of a square root can be cancelled out, as the value for that dimension is also the magnitude of the whole vector, but a 1D intercept occurring on a line is a bit trivial, don't you think?). The way to eliminate the square root is to square both sides. If Tt = Pt it follows that Tt2 = Pt2 and proving this is as simple as referencing the reflexive property after n * n = n * n becomes n = n after dividing both sides by n. We have succeeded in eliminating the square root, but now we must resolve all of the squares we have applied to the sub-equations.

Here is what we have so far: Xt2 + Yt2 + Zt2 = Pt2 will be true for values of t which are the length of time that passes for a straight line intercept (and therefore an optimum intercept, as the fastest way between two points is a straight line. It follows that the fastest, i.e. smallest value of t that is real, greater than or equal to zero, and less than or equal to the remaining burn time is the best intercept, as it is the fastest straight line, the fastest of the fastest, if you will).

To complete this problem, one must internally square each sub-equation. Fortunately the workload is reduced by the fact that the equation is essentially the same for all of them. Later when I mention the calculation of intercepts where the missile runs out of fuel during its flight, the P function will be slightly different, but here it closely resembles the X, Y, and Z functions as well. Because applying the square to the function is so similar, I will only solve for X here, rather than wasting space by copying and pasting the work three more times and changing some letters.

We start with this:

Xt2 = ((1/2) * Tax * t2 + Tvx * t + Tpx)2 

It follows that:

(A + B + C)2 = A * A + A * B + A * C +
               B * A + B * B + B * C +
               C * A + C * B + C * C 

This is similar to the old FOIL method of multiplying binomials (First, Outside, Inside, Last), which is far more common, but it is a broader implementation of the underlying mathematical truth (this one uses trinomials) and it can be simplified for our purposes because we are only interested in squaring trinomials, so we can simplify that:

(A + B + C)2 = A2 + B2 + C2 + 2 * (A * B + A * C + B * C)

Ergo, in standard polynomial form for t with all like terms combined,

Xt2 = (1/4) * Tax2 * t4 +
  Tax * Tvx * t3 +
  (Tax * Tpx + Tvx2) * t2 +
  2 * Tvx * Tpx * t +
  Tpx2

Combining the X, Y, and Z functions is as simple as combining like terms further (the first term, i.e. the t4 term becomes: (1/4) * (Tax2 + Tay2 + Taz2) * t4).

So now, how does this help us solve for t? By subtracting Pt2 from both sides of our original Tt2 = Pt2 (thereby yielding Tt2 - Pt2 = 0), the result is, when everything is fully substituted and reduced in these sub-equations, discover a quartic, i.e. a fourth degree polynomial. To actually perform this last step, subtraction of the Pt2 terms from the Tt2 equation that we completed in the last two paragraphs is as simple as it was to add terms in the last paragraph for X, Y, and Z. Here is the entire resultant equation (I have added some parenthesis to make some things clearer that are not truly necessary for order of operations):

0 = ( (1/4) * (Tax2 + Tay2 + Taz2 - Pam2) ) * t4 +
    ( Tax * Tvx + Tay * Tvy + Taz * Tvz - Pam * Pvm ) * t3 +
    ( Tax * Tpx + Tvx2 + Tay * Tpy + Tvy2 + 
        Taz * Tpz + Tvz2 - Pam * Ppm - Pvm2 ) * t2 +
    ( 2 * ( Tvx * Tpx + Tvy * Tpy + Tvz * Tpz - Pvm * Ppm ) ) * t +
    ( Tpx2 + Tpy2 + Tpz2 - Ppm2 )

This can be solved for t. I did some research and programmed in a rootsOf function, where the coefficients of t are passed in the order of standard polynomial form, and all real roots are returned. To be fair, it isn't a truly "complete" roots function because it simply tells the caller the real roots, but gives no indication of complex roots (they are omitted all together, i.e. if all roots are complex, the function does not return any roots). While this is true, it is more than wishful thinking to say that one has locked onto their target because they can intercept a target at an imaginary time, this is no more valid than an intercept occurring at a negative time for our purposes, but negative numbers are thrown out higher in the call stack within the guidance system, along with other unacceptable values (roots after rbt found with this equation, for instance.)

To find out what our roots are for times after rbt, the only thing that changes is our P function, which then looks like this:

P't = ( Pam * Rbt *+ Pvm ) * t + ( Ppm - (1/2) * Pam * Rbt2 )

which squared looks like this:

P't2 = ( Pam2 * Rbt2 + 2 * Pam * Pvm * Rbt + Pvm2 ) * t2 +
  ( 2 * Pam * Ppm * (Pam * Rbt + Pvm) - Pam2 * Rbt3 - Pam * Pvm * Rbt2 ) * t +
  ( (1/4) * Pam2 * Rbt4 - Pam * Ppm * Rbt2 + Ppm2 )

Please bear in mind that, after subtracting these terms from the T function, the resulting roots are considered valid if they are both real and greater than rbt (roots exactly equal to rbt are also technically valid, but the P function handles these just fine, so there's no need to keep them if they show up when using the P function instead. Please also note that being greater than rbt also means being greater than or equal to zero, as rbt cannot be negative in context).

Lastly, a word about impacts which are glancing blows: Given a quartic polynomial Y = AX4 + BX3 + CX2 + DX + E, it follows that the values for X at which the instantaneous slope is zero, are given by the roots of the cubic polynomial Y' = 4AX3 + 3BX2 + 2CX + D. Please note that the implied coefficient of D is one, and the implied power of X of the D term is zero, which simplifies to D, but also note that if you continue the pattern, the coefficient of E is zero, and this is why the E term disappears altogether.

By using this cubic, we can take our quartic coefficients and calculate closest approaches. The potentially confusing part about this is the fact that they are not necessarily meaningful, as they are commonly "furthest" approaches if a true intercepts or true closest approaches exist at the base of that apex. However, if in fact a closest approach does exist (for instance, if no true intercepts exist, then either the first and last cubic roots are closest approaches, or the middle root is a closest approach, where order is determined by the value of X), then we can determine that it is a valid glancing blow impact if the absolute value of Y' is less than or equal to the radius of tolerance (the radius of the target plus the radius of the missile). Weather an approach is closest or furthest is irrelevant to the guidance system, because the radius of tolerance comparison will eliminate unacceptable values, and whether it is technically a closest or furthest approach is really irrelevant because if the approach is close enough to be an impact, and is also faster than a direct impact, firing at the glancing blow impact is best.

A generic cubic:

A generic quartic:

A cubic with roots at the places where the slope of a quartic is zero:

So, we've solved for all of our roots, now here comes the final step. By substituting those values of t into the X, Y, and Z functions the result is the X, Y, and Z coordinates of the target relative to the shooter (with the shooter's acceleration being completely disregarded, as mentioned long ago). The guidance system can compile an array of unit vectors by taking each of these (X, Y, Z) coordinates and dividing all components and dividing by the respective magnitude of that particular coordinate, again given by M = √(X2 + Y2 + Z2), i.e. the unit vector for coordinate C is given by (X, Y, Z) where X = Cx / Cm Y = Cy / Cm and Z = Cz / Cm.

To use each of these unit vectors, the missile's initial state is set as follows (Acceleration, Velocity, and Position), given a targeting unit vector U and the shooter vectors for velocity and position Sv and Sp (note Pam, Pvm, and Ppm are multiplied via scalar multiplication and that Sv and Sp are added with vector addition):

A = U * Pam
V = U * Pvm + Sv
P = U * Ppm + Sp

Note that during the flight when the missile is updating it's trajectory, it tells the guidance system that it's current velocity and position are the "shooter" parameters, and that it has an initial velocity of zero and an initial position of zero, because adjustments to the flight path must be done at this point with the missile's acceleration alone, i.e. there is no "shooter" to give it an initial velocity or position from where the missile starts out.

Where the Quartic Formula Comes from and What Makes it Significant

In 1540, Lodovico Ferrari discovered the quartic formula. This masterful formula was first published in the 1545 publication entitled Ars Magna, along with the cubic formula, discovered by Gerolamo Cardano, which is so crucial to solving quartic formula that it caused this five year delay. An important question arises—what is a quartic? If there are five defined constants and one variable in an equation, it is quartic if zero equals the first constant times the variable raised to the fourth power, plus the second constant times the variable raised to the third power, plus the third constant raised to the second power, plus the fourth constant times the variable, plus the fifth constant. This is generally written as 0 = AX4 + BX3 + CX2 + DX + E. The quartic formula solves for all four possible values of X that satisfy this condition, where A != 0. There are always technically four solutions, but, before removing repeated solutions, zero, two, or four of them will be real numbers, and the remaining solutions for X will be complex numbers—numbers containing a real part plus an imaginary part, the latter of which is proceeded by the coefficient i. These four solutions are extremely useful in a wide variety of real world situations, because having to square a squared term to come to a conclusion, and thereby taking that term to the fourth power, is surprisingly common in practice.

This is particularly insightful to a guided missile that takes the acceleration, velocity, and position of itself and its target into account. After substituting, distributing, combining like terms, and methodically moving everything to one side of the equation, hence leaving zero on the other side, a quartic equation results, due to the fact that acceleration changes position in respect to time appears as pt = (1/2)at2, and in the process of coming to the quartic representation, this becomes squared and appears as pt2 = (1/4)a2t4, where we wish to solve for t, rather than X, the A term unfolds before our eyes. In a similar fashion, a guided missile that only takes velocity and position into account will produce results through use of a quadratic equation, generally written as 0 = AX2 + BX + C, because in the process of aligning everything in standard polynomial form t will be squared, but because there are no t squared terms to be squared, the highest power of t in the resulting form will be two, hence AX2. On a yet simpler level, a guided missile that only takes position into account will produce even a simpler solution.

These conclusions suggest that a missile might be able to navigate with knowledge of position, velocity, acceleration, and change in acceleration, in the same sense as acceleration represents a change in velocity, though there is no good descriptor of a change in acceleration. Surprisingly enough, this is, however, incorrect. Another question arises—what makes acceleration so special? The answer lies in the power of four that results in a quartic itself. A quintic, or an equation of any nonzero coefficient of a higher degree polynomial, does not have a formula that, given the coefficients, produces all possible values for the variable in question, namely time in the case of the guided missile. It is common knowledge in the mathematical community that "The Abel–Ruffini theorem (also known as Abel's impossibility theorem) states that there is no general solution in radicals to polynomial equations of degree five or higher" (Wikipedia). There are methods for homing in on really good estimations of all of the roots of such a polynomial, but there is no way to use the coefficients in an elegant equality statement to find all of its roots as the quadratic, cubic, and quartic formulas provide for polynomials of the second, third, and fourth degree.

Points of Interest

I am aware that as fuel is exhausted the mass of the projectile will decrease, and, applying F = MA (Newton's second law of motion) to the problem, we find that we have produced change in acceleration. Unfortunately I could generate a higher degree polynomial, but I could not solve it analytically for time, as the quartic formula is the highest degree polynomial that can be solved in this way (the Wikipedia reference in the above section details this at length), and because this motion is relatively speaking quite small, and can be reasonably ignored. This guidance system is for situations where there are no forces such as air resistance or gravity, but again the inclusion of such forces would cause relatively small errors over time, as the missile would constantly update its trajectory at every time-step to account for these small adjustments. I am still working on broader solutions to the problem, but my findings are large enough at this point (and functional enough too) that I thought I might as well share.

History

First I made a rough routine for estimation of constant velocities in 2D in a matter of hours, later I developed an analytical 2D solution using trigonometry (primarily the law of cosines, solving for some unknowns with the quadratic equation). I was able to convert this into a 3D solution for constant velocities, but hit a wall here. After changing my thinking and moving towards algebra and calculus, and away from geometry and trigonometry, I have a 3D solution for constant acceleration, which is easily adaptable for any number of dimensions, as the quartic coefficients are quite simply the sum of the results of a particular small sub-equation for every dimension. Want a 6D solution? how about an 11D solution? one can easily set that up in 5 minutes using this method of solving for time.

• 22 June, 2007 -- Original version posted
• 27 June, 2007 -- Article and downloads updated
• 03 July, 2007 -- Article completed and simple UI download added
• 10 July, 2007 -- Download Simple_Guidance_System_UI.zip updated
• 20 July, 2007 -- Download Simple_Guidance_System_UI.zip updated
• 7 August, 2007 -- Download Simple_Guidence_System_UI.zip updated
• 9 August, 2007 -- Download Simple_Guidence_System_UI.zip updated
• 15 August, 2007 -- Bug fixed in the download Guided_Missile_Simulator.zip that prevented some valid user input from passing through as valid