Simulator Overview

This simulator is designed to allow us to explore strategies and tactics for swarm-on-swarm combat, particularly in the area of missile and UAV defense, and with a bias toward small smart munitions. Our aspirational goal is to develop a science of swarm-on-swarm combat. In phase 1 of the project we have implemented a simple aerodynamic model that provides rough estimates of the capabilities of different weapons. Proportional navigation and augmented proportional navigation guide interceptors to threats. Threats have the ability to evade. Threats and interceptors have perfect position knowledge at a configurable sensor update rate. At this stage, interceptor launch and submunition dispense are hand coded.

Moving forward we will add sensor range and resolution limits, and a realistic communication model. Future versions will explore optimal control and machine learning approaches to launch sequencing, target assignment, trajectory generation, and control.

Introduction

The simulator performs a multi-agent simulation between two types of agents: interceptors and threats. The threats will target the static asset, located at the origin of the coordinate system, and the interceptors will defend the asset from the incoming threats.

There are two types of interceptors:

• Carrier interceptors: interceptors that carry and dispense other interceptors (e.g., Hydra-70 rockets)
• Missile interceptors: interceptors that pursue threats (e.g., micromissiles)

There are also two types of threats:

• Fixed-wing threats: Pursue their targets using proportional navigation (PN)
• Rotary-wing threats: Pursue their targets using direct linear guidance

Simulator Physics

Agent System Model

Each agent is modeled as a point mass, i.e., a 3-DOF body ignoring any rotations. It also has instantaneous acceleration in all directions, subject to constraints, because we do not model any actuator or rotational dynamics. Finally, we abstract away the aerodynamics of the agents, so we do not model the angle of attack or stall.

As a point mass, each agent is represented by a six-dimensional state vector consisting of the agent's three-dimensional position and three-dimensional velocity. The input to the system is a three-dimensional acceleration vector.

The state vector is given by:

$$\overrightarrow{x}(t)=\left[\begin{array}{c}\overrightarrow{p}(t)\\ \overrightarrow{v}(t)\end{array}\right]\in {\mathbb{R}}^6,$$

where $\overrightarrow{p}(t)\in {\mathbb{R}}^3$ denotes the agent's position and $\overrightarrow{v}(t)\in {\mathbb{R}}^3$ denotes the agent's velocity in the Cartesian coordinates.

The input vector is given by:

$$\overrightarrow{u}(t)=\overrightarrow{a}(t)\in {\mathbb{R}}^3,$$

where $\overrightarrow{a}(t)\in {\mathbb{R}}^3$ denotes the agent's acceleration.

The nonlinear state evolution equation is given by:

$$\frac{d}{dt}\overrightarrow{x}(t)=\left[\begin{array}{c}\overrightarrow{v}(t)\\ \overrightarrow{a}(t)-\overrightarrow{g}-(\frac{F_D(\overrightarrow{v}(t))}{m}+\frac{\Vert {\overrightarrow{a}}_{\perp }(t)+{\overrightarrow{g}}_{\perp }\Vert }{(L/D)})\frac{\overrightarrow{v}(t)}{\parallel \overrightarrow{v}(t)\parallel }\end{array}\right],$$

where $\overrightarrow{g}=\left[\begin{array}{c}0\\ 0\\ g\end{array}\right]$ represents the acceleration due to gravity, $\frac{F_D(\overrightarrow{v}(t))}{m}$ represents the deceleration along the agent's velocity vector due to air drag, and $\frac{\Vert {\overrightarrow{a}}_{\perp }(t)+{\overrightarrow{g}}_{\perp }\Vert }{(L/D)}$ represents the deleceration along the agent's velocity vector due to lift-induced drag.

Any acceleration normal to the agent's velocity vector, including the components of the acceleration vector ${\overrightarrow{a}}_{\perp }(t)$ and gravity vector ${\overrightarrow{g}}_{\perp }$ that are normal to the velocity vector, will induce some lift-induced drag.

$${\overrightarrow{a}}_{\perp }(t)=\overrightarrow{a}(t)-{\text{proj}}_{\overrightarrow{v}(t)}(\overrightarrow{a}(t))$$$${\overrightarrow{g}}_{\perp }=\overrightarrow{g}-{\text{proj}}_{\overrightarrow{v}(t)}(\overrightarrow{g})$$

Agent Acceleration

The agent acceleration is given by:

$$\frac{d}{dt}\overrightarrow{v}(t)=\overrightarrow{a}(t)-\overrightarrow{g}-(\frac{F_D(\overrightarrow{v}(t))}{m}+\frac{\Vert {\overrightarrow{a}}_{\perp }(t)+{\overrightarrow{g}}_{\perp }\Vert }{(L/D)})\frac{\overrightarrow{v}(t)}{\parallel \overrightarrow{v}(t)\parallel }$$

Unlike interceptors, threats are not subject to drag or gravity.

The air drag is given by:

$$F_D(\overrightarrow{v}(t))=\frac{1}{2}\rho C_{D}A\parallel \overrightarrow{v}(t){\parallel }^2,$$

where $\rho$ is the air density that decays exponentially with altitude: $\rho =1.204\frac{\text{kg}}{{\text{m}}^3}\cdot e^{-\frac{\text{altitude}}{10.4\text{ km}}}$, $C_D$ is the airframe's coefficient of drag, and $A$ is the cross-sectional area. For all angles of attack, we specify a constant $(L/D)$ ratio.

We do impose some constraints on the acceleration:

• Interceptors can only accelerate normal to their velocity (no thrust during the midcourse phase), i.e., $\overrightarrow{a}(t)\cdot \overrightarrow{v}(t)=0$. Therefore, $\overrightarrow{a}(t)={\overrightarrow{a}}_{\perp }(t)$ for interceptors.
• Threats may have some forward acceleration, which is bounded by the maximum forward acceleration specified for each threat type.
• The normal acceleration is constrained by the maximum number of g's that the agent's airframe can pull:$$\parallel {\overrightarrow{a}}_{\perp }(t)\parallel \le {\left(\frac{\parallel \overrightarrow{v}(t)\parallel }{v_{\text{ref}}}\right)}^2{a}_{\text{ref}}$$$a_{\text{ref}}$ denotes the maximum normal acceleration that the airframe can pull at the reference speed $v_{\text{ref}}$.

Simulator Behaviors

Sensing

Currently, all agents are equipped with an ideal sensor, one that can peek through the fog of war with no noise and no delay. Sensing is performed within the agent's frame of reference using spherical coordinates, so each sensor output $\overrightarrow{y}$ is a nine-dimensional vector.

$$\overrightarrow{y}=\left[\begin{array}{c}{\overrightarrow{y}}_p\\ {\overrightarrow{y}}_v\\ {\overrightarrow{y}}_a\end{array}\right]\in {\mathbb{R}}^9,$$

where ${\overrightarrow{y}}_p\in {\mathbb{R}}^3$ denotes the three-dimensional position difference between the agent and its sensing target, ${\overrightarrow{y}}_v$ denotes the three-dimensional velocity difference between the agent and its sensing target, and ${\overrightarrow{y}}_a$ denotes the three-dimensional acceleration of the sensing target. ${\overrightarrow{y}}_p$ and ${\overrightarrow{y}}_v$ are both given in spherical coordinates in the agent's frame of reference while ${\overrightarrow{y}}_a$ is in Cartesian coordinates.

Interceptors are constrained in their sensor update frequency, which is configurable for each interceptor type. As a result, interceptors can change their actuation input at a rate faster than the sensor update frequency. The simulator currently uses a naive guidance filter simply performs a zero-hold interpolation on the latest sensor output and applies the latest acceleration to a model of the sensing target until the next sensor output arrives and resets the model's position, velocity, and acceleration. In other words, for $nT\le t<(n+1)T$, where $T$ is the sensor update period, the simple target model is as follows:

$$\frac{d}{dt}\left[\begin{array}{c}\overrightarrow{p}(t)\\ \overrightarrow{v}(t)\end{array}\right]=\left[\begin{array}{c}\overrightarrow{v}(t)\\ \overrightarrow{a}(t){|}_{t=nT}\end{array}\right],$$

where the initial conditions $\overrightarrow{p}(t){|}_{t=nT}$ and $\overrightarrow{v}(t){|}_{t=nT}$ are set by the latest sensor output.

Threats are assumed to be omniscient, so they have no frequency constraint on their sensor output and know the positions, velocities, and accelerations of all interceptors at all times.

Guidance & Navigation

Proportional Navigation

Using the fact that constant bearing decreasing range (CBDR) leads to a collision, we apply an acceleration normal to the velocity vector to correct for any bearing drift. In the simulator, proportional navigation follows the simple control law:

$${\overrightarrow{a}}_{\perp }=K\dot{\overrightarrow{\lambda }}v,$$

where $K$ is the navigation gain, $\dot{\overrightarrow{\lambda }}$ is the rate of change of the bearing, and $v$ is the closing velocity. For interceptors, we choose $K=3$.

Proportional navigation is effective for non-accelerating targets and guarantees a collision. However, simply using true proportional navigation as a guidance law leads to some undesired behavior when the rate of change of the bearing $\dot{\overrightarrow{\lambda }}$ is near zero.

1. Increasing range: The closing velocity may be negative, i.e., the distance between the agent and its target may actually be increasing. In this case, the agent should apply a maximum normal acceleration in any direction to turn around, but since $\parallel \dot{\overrightarrow{\lambda }}\parallel \approx 0$, the normal acceleration is minimal. To overcome this issue, the navigation gain is increased significantly if the closing velocity is negative.

   // Handle the case where the closing velocity is negative.
   if (closingVelocity < 0) {
     navigationGain = Mathf.Max(1f, Mathf.Abs(closingVelocity) * 100f);
   }

2. Spiral behavior: If the target is at a near-constant ${90}^{\circ }$ bearing from the agent, the agent may end up in a spiral encircling the target because $\overrightarrow{\lambda }$ is roughly constant and so $\parallel \dot{\overrightarrow{\lambda }}\parallel \approx 0$. To overcome this limitation, the agent will apply a higher navigation gain if the target bearing is within $\pm {10}^{\circ }$ of ${90}^{\circ }$.

   // Handle the case of the spiral behavior if the target is at a bearing of 90 degrees +- 10 degrees.
   if (Mathf.Abs(Mathf.Abs(sensorOutput.position.azimuth) - Mathf.PI / 2) < 0.2f ||
       Mathf.Abs(Mathf.Abs(sensorOutput.position.elevation) - Mathf.PI / 2) < 0.2f) {
     // Clamp the line-of-sight rate at 0.2 rad/s.
     float minLosRate = 0.2f;
     losRateAz = Mathf.Sign(losRateAz) * Mathf.Max(Mathf.Abs(losRateAz), minLosRate);
     losRateEl = Mathf.Sign(losRateEl) * Mathf.Max(Mathf.Abs(losRateEl), minLosRate);
     navigationGain = Mathf.Abs(closingVelocity) * 100f;
   }

Augmented Proportional Navigation

Augmented proportional navigation adds a feedthrough term proportional to the agent’s acceleration:

$${\overrightarrow{a}}_{\perp }=K(\dot{\overrightarrow{\lambda }}v+\frac{1}{2}{\overrightarrow{a}}_{\text{target}}),$$

where ${\overrightarrow{a}}_{\text{target}}$ is the target’s acceleration that is normal to the agent's velocity vector.

Augmented proportional navigation is equivalent to true proportional navigation if the target is not accelerating.

Ground Avoidance

Gravity only acts on interceptors as the simulator assumes that the threats are able to compensate for gravity. Currently, interceptors do not consider gravity when determining their acceleration input for the next simulation step. As a result, gravity acts as a disturbance to each interceptor's dynamics system and may cause the interceptor to collide with the ground if not accounted for.

Threats may also collide with the ground, especially after having performed an evasion maneuver from pursuing interceptors. The simulator implements a basic ground avoidance algorithm for the threats: If the threat's vertical speed will cause the threat to collide with the ground before the threat will hit the asset, the threat will adjust its vertical velocity to be a linear combination of navigating towards the asset and pulling up away from the ground.

// Calculate the time to ground.
float altitude = transformPosition.y;
float sinkRate = -transformVelocity.y;
float timeToGround = altitude / sinkRate;

// Calculate the time to target.
float distanceToTarget = sensorOutput.position.range;
float timeToTarget = distanceToTarget / transformSpeed;

float groundProximityThreshold = Mathf.Abs(transformVelocity.y) * 5f;
if (sinkRate > 0 && timeToGround < timeToTarget) {
  // Evade upwards normal to the velocity.
  Vector3 upwardsDirection = Vector3.Cross(transformForward, transformRight);

  // For the y-component, interpolate between the calculated acceleration input and the upward acceleration.
  float blendFactor = 1 - (altitude / groundProximityThreshold);
  accelerationInput.y = Vector3.Lerp(accelerationInput, upwardsDirection * CalculateMaxNormalAcceleration(), blendFactor).y;
}

Interceptor Assignment

Threat-Based Assignment

When submunitions are dispensed (e.g., micromissiles from Hydra-70s), they are assigned a threat to intercept. The assignment algorithm prefers assigning an interceptor to each threat before doubling up on previously assigned threats.

Given the list of threats, the simulator first sorts the threats as follows:

1. Sorting (descending) by the number of already assigned interceptors
2. Sorting (ascending) by threat value, where the threat value of a threat is given by:$$V_{threat}=\frac{\parallel \overrightarrow{v}(t)\parallel }{\parallel \overrightarrow{p}(t)\parallel },$$where $\parallel v_t\parallel$ is the threat's speed and $\overrightarrow{p}(t)$ is the threat's position vector, assuming that the asset is at the origin. After sorting the threats, we simply assign interceptors down the list.

Note that this algorithm may not be optimal but is a good starting point.

Intercept Evasion

When interceptors get too close to their intended target, the threat performs an evasive maneuver to expend the interceptor's speed and remaining energy. During the evasive maneuver, the threat performs the following:

1. The threat accelerates to its maximum speed.
2. The threat turns away from the incoming interceptor at its maximum normal acceleration and tries to align its velocity vector to be normal to the interceptor's velocity vector. Since the threat applies a normal acceleration, the interceptor must turn too and thus sacrifice speed due to the lift-induced drag.

// Evade the interceptor by aligning the threat's velocity vector to be normal to the interceptor's velocity vector.
Vector3 normalVelocity = Vector3.ProjectOnPlane(transformVelocity, interceptorVelocity);
Vector3 normalAccelerationDirection = Vector3.ProjectOnPlane(normalVelocity, transformVelocity).normalized;

// Turn away from the interceptor.
Vector3 relativePosition = interceptorPosition - transformPosition;
if (Vector3.Dot(relativePosition, normalAccelerationDirection) > 0) {
  normalAccelerationDirection *= -1;
}

If the threat is too close to the ground, however, it must ensure that it does not collide with the ground. Therefore, as it approaches the ground, the threat instead performs a linear combination of:

1. turning to evade the interceptor, as described above, and
2. turning parallel to the ground.

// Avoid evading straight down when near the ground.
float altitude = transformPosition.y;
float sinkRate = -transformVelocity.y;
float groundProximityThreshold = Mathf.Abs(transformVelocity.y) * 5f;
if (sinkRate > 0 && altitude < groundProximityThreshold) {
  // Determine evasion direction based on the bearing to the interceptor.
  Vector3 toInterceptor = interceptorPosition - transformPosition;
  Vector3 rightDirection = Vector3.Cross(Vector3.up, transform.forward);
  float angle = Vector3.SignedAngle(transform.forward, toInterceptor, Vector3.up);

  // Choose the direction that turns away from the interceptor.
  Vector3 bestHorizontalDirection = angle > 0 ? -rightDirection : rightDirection;

  // Interpolate between horizontal evasion and upward ground avoidance.
  float blendFactor = 1 - (altitude / groundProximityThreshold);
  normalAccelerationDirection = Vector3.Lerp(normalAccelerationDirection, bestHorizontalDirection + transform.up * 5f, blendFactor).normalized;
  }