Impact distribution¶
$$ \newcommand{\eg}{{\it e.g.}} \newcommand{\ie}{{\it i.e.}} \newcommand{\argmin}{\operatornamewithlimits{argmin}} \newcommand{\mc}{\mathcal} \newcommand{\mb}{\mathbb} \newcommand{\mf}{\mathbf} \newcommand{\minimize}{{\text{minimize}}} \newcommand{\diag}{{\text{diag}}} \newcommand{\cond}{{\text{cond}}} \newcommand{\rank}{{\text{rank }}} \newcommand{\range}{{\mathcal{R}}} \newcommand{\null}{{\mathcal{N}}} \newcommand{\tr}{{\text{trace}}} \newcommand{\dom}{{\text{dom}}} \newcommand{\dist}{{\text{dist}}} \newcommand{\R}{\mathbf{R}} \newcommand{\SM}{\mathbf{S}} \newcommand{\ball}{\mathcal{B}} \newcommand{\bmat}[1]{\begin{bmatrix}#1\end{bmatrix}} \newcommand{\loss}{\ell} \newcommand{\eloss}{\mc{L}} \newcommand{\abs}[1]{| #1 |} \newcommand{\norm}[1]{\| #1 \|} \newcommand{\tp}{T} $$In this problem, the objective is to simulate the trajectory of a free-flying rocket and estimate the dispersion of the impact point due to wind effects.
Let us consider the equations of motion for a free-flying rocket.
The dynamics of the state variables, which include the velocity (speed), $V$, the flight path angle in the vertical plane, $\gamma$, and the heading angle in the horizontal plane, $\psi$, are given as follows.
The vertical flight path angle, $\gamma$ is defined as 0 degrees when the rocket is flying horizontally, increasing as the rocket ascends, with a value of 90 degrees representing vertical ascent. The horizontal heading angle is defined as 0 degrees when flying to north, increasing as the rocket turns right, reaching 90 degrees when flying to east.
\begin{align*} m\dot{V} &= D_x - mg\sin\gamma \\ mV\dot{\gamma} &= -D_z -mg\cos\gamma \\ mV\dot{\psi} &= D_y \\ \end{align*}
Here, $m$ represents the mass of the rocket, and the drag force $D = (D_x, D_y, D_z) \in \mathbb{R}^3$ consists of three components: 1) $D_x$, the drag component along the direction of the velocity vector, 2) $D_y$, the horizontal drag component perpendicular to the velocity vector and to the right, and 3) $D_z$, the vertical drag component perpendicular to the velocity vector and acting downward. In general, a free-flying rocket is assumed to be statically stable, meaning the angle of attack remains 0 degrees at all times, allowing us to assume that no lift is generated.
The position of the rocket, relative to the launch point, is expressed as $p = (p_n, p_e, h) \in \mathbb{R}^3$, where $p_n$ represents the position in the northward ($N$) direction, $p_e$ represents the displacement in the eastward ($E$) direction, and $h$ represents the altitude ($H$). Note that the altitude ($H$) is in the opposite direction of the downward ($D$) direction, where gravity acts.
The rate of the position variables $p_n$, $p_e$, and $h$ is described by the following dynamics, where the right-hand sides can be interpreted as the $NEH$ components of the velocity vector, $v_n$ (northward velocity), $v_e$ (eastward velocity), and $\dot{h} = -v_d$(velocity in the altitude direction), expressed in terms of the state variables $V$, $\gamma$, and $\psi$.
\begin{align*} \dot{p}_n &= v_n = V\cos\gamma\cos\psi \\ \dot{p}_e &= v_e = V\cos\gamma\sin\psi \\ \dot{h} &= - v_d = V\sin\gamma \end{align*}
The dynamics of the rocket’s trajectory can be expressed using the six state variables described above. With these six variables, the full equations of motion can be formulated. Once the drag components $D = (D_x, D_y, D_z)$, acting in the direction of motion and perpendicular directions, are provided at each moment, the time derivatives of the six state variables can be computed. These time derivatives can then be integrated over time to simulate the rocket’s trajectory.
Now, let us examine how the drag force is determined. The magnitude of the drag force is proportional to the dynamic pressure, which is modeled as being proportional to the atmospheric density (a function of altitude) and the square of the rocket’s relative velocity with respect to the wind. The direction of the drag force is opposite to the relative velocity to the wind. The parameters $C_d$ and $S$ represent the drag coefficient and the cross-sectional area of the rocket, respectively, and are assumed to be constant.
Combining these factors, the drag components in the North-East-Down (NED) coordinate system, $D_{ned} = (D_n, D_e, D_d)$, can be computed as follows:
\begin{align*} D_{ned} = \bmat{D_n \\ D_e \\ D_d} &= \overbrace{\overbrace{\frac{1}{2}\rho \|v_\text{rel} \|^2}^{\text{Dynamic pressure}}SC_D}^{\text{Magnitude of drag}} \overbrace{\left(-\frac{v_\text{rel}}{\|v_\text{rel} \|}\right)}^{\text{Direction of drag}} \\ &= -\frac{1}{2}\rho SC_D \|v_\text{rel} \| v_\text{rel} \end{align*}
The atmospheric density $\rho$ can be modeled as a function of altitude $h$ as follows, where the density is in $\text{kg/m}^3$, and the altitude is in $\text{m}$:
$$ \rho(h) = \rho_0 \exp\left(-\frac{h}{H_s}\right) $$
where:
- $\rho_0$ is the atmospheric density at sea level (approximately $1.225 \text{kg/m}^3$),
- $H_s$ is the scale height, which is typically around $8,500 \text{m}$ in the Earth’s atmosphere,
- $h$ is the altitude in meters.
This exponential model accounts for the decrease in atmospheric density with increasing altitude.
The relative velocity of the rocket with respect to the wind, $v_{\text{rel}} \in \mathbb{R}^3$, can be expressed in the North-East-Down (NED) coordinate system as follows:
\begin{align*} v_\text{rel} &= v_{ned}- v_{w} \\ &= \bmat{V\cos\gamma\cos\psi \\ V\cos\gamma\sin\psi \\ -V\sin\gamma} - \bmat{v_{w,n} \\ v_{w,e} \\ v_{w,d}} \end{align*}
The drag components in the NED coordinate system, $D_{ned} = (D_n, D_e, D_d)$, can be transformed into the XYZ coordinate system which aligns with the velocity vector, yielding the XYZ drag components $D_{xyz} = (D_x, D_y, D_z)$. This transformation is performed using a rotation matrix that converts the NED frame into the XYZ frame based on the flight path angle $\gamma$ and heading angle $\psi$.
\begin{align*} D_{xyz} = \bmat{D_x \\ D_y \\ D_z} = \bmat{\cos\gamma\cos\psi & \cos\gamma\sin\psi & -\sin\gamma \\ -\sin\psi & \cos\psi & 0 \\ \sin\gamma\cos\psi & \sin\gamma\sin\psi & \cos\gamma } D_{ned} \end{align*}
The simulation parameters are defined as follows.
- $m$: Mass of the rocket ($=40\,\text{kg}$)
- $d$: Diameter of the rocket ($=0.16\,\text{m}$)
- $S$: Cross-sectional area of the rocket ($=\frac{\pi d^2}{4}$)
- $C_d$: Drag coefficient of the rocket ($=0.2$)
- $g$: Gravitational acceleration ($=9.8\,\text{m/s}^2$)
These constants are used to calculate drag forces and the overall dynamics of the rocket throughout its trajectory.
import numpy as np
import matplotlib.pyplot as plt
# parameters
m = 40 # mass
S = np.pi*0.08**2 # cross-section area
Cd = 0.2 # drag coefficient
g = 9.8 # gravitational acceleration
# unit conversion
Rad2Deg = 180/np.pi # conversion: radian to degree
Deg2Rad = 1/Rad2Deg # conversion: degree to radian
(Problem 1) In this problem, we aim to compute and visualize the 3D trajectory of a rocket under the assumption that there is no wind ($v_{w} = 0$).
The trajectory will be integrated using numerical methods. One can perform the numerical integration using methods such as the Adams-Bashforth technique or alternatively by leveraging functions provided in libraries like scipy.
The goal is to numerically integrate the equations of motion and plot the 3D trajectory of the rocket based on the dynamics described earlier. The absence of wind simplifies the drag calculations, allowing us to focus on the effects of the rocket’s motion and aerodynamic drag in a wind-free environment.
(Problem 1a) Let the state vector z of the system be defined as:
$$ z = \begin{bmatrix} V & \gamma & \psi & p_n & p_e & h \end{bmatrix}^T $$and use the equations of motion without considering the wind, to define a function state_derivative that calculates the time derivatives $\dot{z}_t$ based on the current state variable $z_t$.
def state_derivative(z):
V, gamma, psi, pn, pe, h = z
##################
# your code here #
# your code here #
##################
return np.array([V_dot, gamma_dot, psi_dot, pn_dot, pe_dot, h_dot])
(Problem 1b) Using the model described above, perform numerical integration of the dynamics of a free-flying rocket launched to east with an initial flight path angle of 40 degrees and an initial velocity of 1 km/s.
Compute and plot the rocket’s altitude-range trajectory plot (altitude in the vertical axis and range in the horizontal axis with range implying the horizontal distance traveled in the launch direction). The numerical integration should be performed using the provided parameters, starting from launch and continuing until the rocket impacts the ground.
# integration parameters
tf = 200 # final time
dt = 0.1 # step size
t = np.arange(0, tf, dt) # time (tentative)
N = len(t) # number of time steps
# initial conditions
V = 1000.0 # velocity
gamma = 40*Deg2Rad # flight path angle
psi = 90*Deg2Rad # heading
pn, pe, h = np.zeros(3) # position
init = [V, gamma, psi, pn, pe, h]
##################
# your code here #
# your code here #
##################
(Problem 1c) Plot the state variables over time using the provided model. Additionally, print out the rocket’s velocity, flight path angle, and flight range (distance traveled) at the moment when it impacts the ground.
##################
# your code here #
# your code here #
##################
Velocity: 395.69 m/s Flight path angle: -58.16 deg Flight range: 33.77 km
(Problem 2) Perform a Monte Carlo simulation to compute the trajectory of a rocket considering the effects of wind. The following assumptions are made:
Wind is assumed to blow only horizontally, meaning there is no vertical component in the $h$- or $d$-direction. In the NED coordinate system, the wind velocity is expressed as $v_w = (v_{w,n}, v_{w,e}, 0)$.
- At the initial time step $k = 0$, the wind velocity components $v_{w,n}$ and $v_{w,e}$ follow a normal distribution with a mean of $0$ and a standard deviation of $5$. The vertical wind velocity component $v_{w,d} = 0$ is set to $0$ for all simulations.
- A random walk is applied to simulate wind speed changes at each time step $k$. The wind velocities are updated as follows, where $r_n$ and $r_e$ are random variables generated independently at each step, both following a normal distribution with a mean of $0$ and a standard deviation of $1$.
(Problem 2a)
Using the same definition of state variables as in the previous problem, define the function state_derivative that computes the time derivatives of the state variables $\dot{z}_t$ considering the wind effects. The function should take the current state $z_t$ and wind velocity $v_w$ as inputs. You may reuse most of the state_derivative function defined in Problem 1, but clearly indicate the modified parts to account for the wind.
def state_derivative(z,v_w):
V, gamma, psi, pn, pe, h = z
##################
# your code here #
# your code here #
##################
return np.array([V_dot, gamma_dot, psi_dot, pn_dot, pe_dot, h_dot])
(Problem 2b) Perform a Monte Carlo simulation with 2000 iterations, each accounting for random wind effects, to simulate the rocket’s trajectory. Overlap all 2000 trajectories on a single plot.
# integration parameters
tf = 200 # final time
dt = 0.1 # step size
t = np.arange(0, tf, dt) # time (tentative)
N = len(t) # number of time steps
# initial conditions
V = 1000.0 # velocity
gamma = 40*Deg2Rad # flight path angle
psi = 90*Deg2Rad # heading
pn, pe, h = np.zeros(3) # position
init = [V, gamma, psi, pn, pe, h]
# Monte-Carlo simulation parameter
n_mc = 2000 # number of Monte-Carlo runs
##################
# your code here #
# your code here #
##################
100 200 300 400 500 600 700 800 900 1000 1100 1200 1300 1400 1500 1600 1700 1800 1900 2000
(Problem 2c)
Using the results of the 2000 Monte Carlo simulations, plot the impact points on the North-East (NE) plane using the scatter function from the matplotlib.pyplot module. Additionally, plot the mean impact point on the same plane.
The accuracy of a rocket’s impact is typically expressed using the Circular Error Probable (CEP). The CEP is defined as the radius of a circle, centered at the mean impact point, within which 50% of the impacts will occur. For example, if a rocket has a CEP of 200 meters, approximately 50 out of 100 launches will land within a 200-meter radius from the mean impact point.
Based on the Monte-Carlo simulation results, estimate the CEP for this rocket.
##################
# your code here #
# your code here #
##################
CEP: 0.70 km
