What is the Circular Error Probable (CEP) for the upcoming drop test of the Electromechanical Recovery System (ERS) under a range of possible scenarios?
Loosely speaking, CEP is the mean distance of a sample of points from their mean impact point. That is, a circle around the average point that contains 50% of the sample points. Here, we are defining CEP as the square root of the mean square error, i.e. the sum of range and azimuth variances and their covariance (ignoring bias for now).
The ERS will be dropped out of a helicopter somewhere near the McMinnville Airport and for safety and liability reasons, we must have an expectation of the area it may land in. The primary variables we are considering is a range of possible helicopter trajectories (in terms of velocity and heading, with positions randomized around a certain mean) as well as a range of likely wind conditions (also as velocities and headings). GPS error is accounted for in the final calculation.
There are four operational modes (that is, three failure modes and a successful mission) to consider. There are also two helicopter airspeeds, 15.5 and 30 knots, to consider, thus we have eight cases of simulations in total.
For our operational modes, we consider the case that the mission is successful. We consider the two cases that only the drogue parachute fails and that the ERS system fails completely to be identical, given the distance we drop from. We consider the case that the drogue deploys correctly but remains attached while the main chute does not deploy. We consider the case that the main chute somehow deploys when the drogue is supposed to, which we expect to provide our worst case in terms of CEP spread.
In all scenarios, we consider the wind in a normal distribution about one direction with a standard deviation of 10 degrees. Thus, all simulations are contained within a 60 degree cone, and a large majority are within a 40 degree cone. The helicopter is considered to be going against that mean wind direction at the specified airspeed with some noise injected for realism. On the plots, the wind blows downwards and the helicopter flies upwards. It should also be noted that we may wish to discard any simulations with wind speeds above 5.15 m/s (that is approximately 10 knots) in order to lower CEP, but we will not do so here.
Since the simulations are conducted relative to an arbitrary direction, similar results should hold in any absolute direction by symmetry. The utility of this approach is that on the day of the test, given the current wind velocity, we should have an idea of how where the system is likely to drift to relative to its release point. This could allow us to select our release point to ensure that the possible impact zone is bounded, and in the case of an irregularly shaped field to drop in, we would have a better sense of which wind speeds are acceptable in any given direction.
Another Notebook contains a trajectory simulator and our documentation of it. We have reason to believe our model is approximately representative of the actual system and environment, assuming wind may be held constant for up to several minutes. The trajectory is available both as a numerical integration of a system of ordinary differential equations and as a linearization of that integration. In either case, the interaction of the drogue parachute with air must be integrated, but the phases prior to that and after that may be considered as steady-states without loss of significant precision. Hence, we use the numerical integration of our trajectory model in this Notebook which is much slower to compute, but much smoother.
Since CEP relies on observational evidence to calculate, we selected a Monte Carlo simulation of our trajectory model to generate data for statistical analysis. For each of our eight scenarios, we randomly sample 1,250 trajectories (with random variables considered independently) and we collect every final impact point, along with their initial conditions.
In [1]:
%run drop_test_simulator.ipynb
import pandas as pd
from bokeh.io import output_file, output_notebook, show
from bokeh.models import (
GMapPlot, GMapOptions, ColumnDataSource, Circle, LogColorMapper, BasicTicker, ColorBar,
DataRange1d, PanTool, WheelZoomTool, BoxSelectTool, Range1d
)
from bokeh.models.mappers import ColorMapper, LinearColorMapper
from bokeh.palettes import Viridis5
In [2]:
def convert_to_lat(y):
lat_in_ft = 364618.1186706298
return y/lat_in_ft
def convert_to_long(x):
long_in_ft = 257742.06079152823
return x/long_in_ft
Note that this code has been written with a lot of generality to make it easier to specify later simulations of our various scenarios completely with single lines of code.
In [3]:
heli_direction = 230
# note that if dt=0.25 the horizontal distance may be off by about 8 m (compared to dt=0.001)
def mc_init(num_sims=500, dt=0.25, mode='', airspeed=15.43334):
all_coords = []
v_wind_rndm = []
v_plane_rndm = []
init_pos_offset = []
deg_wind_mean, deg_wind_std_dev = (50, 10) # 20-40-60 degree cone (with probability 68-95-99.7)
for i in range(num_sims):
v_wind_rndm.append([np.random.rayleigh(2.47344), # rayleigh distro good for wind magnitude
np.random.normal(deg_wind_mean, deg_wind_std_dev)])
v_plane_rndm.append(np.random.normal(
airspeed + v_wind_rndm[-1][0] * np.cos(np.radians(v_wind_rndm[-1][1] - heli_direction)),
1))
init_pos_offset.append([np.random.normal(0, 30 * 0.3048), # converting feet to meters
np.random.normal(0, 30 * 0.3048),
np.random.normal(0, 15 * 0.3048)]) # account for the imperfections of reality
print('random wind and helicopter vectors obtained\n')
for i, vw in enumerate(v_wind_rndm):
all_coords.append(trajectory(v_plane=v_plane_rndm[i],
deg_plane=heli_direction,
v_wind=vw[0],
deg_wind=vw[1],
x_0_offset=init_pos_offset[i],
mode=mode,
dt=dt)[0])
if i % 250 == 0: print('iterations:', i)
print('done simulations!\n')
return (all_coords, v_wind_rndm, v_plane_rndm)
In [4]:
def polar_coord(df):
# magnitude of landing coordinate
df['norm'] = df.loc[:,['x','y']].apply(np.linalg.norm, axis=1, raw=True)
# azimuth angle of landing (note that this is on the unit circle)
df['theta'] = df.apply(lambda pos: np.arctan2(pos[1], pos[0]), axis=1, raw=True)
def lat_long(df):
df['latitude'] = df.apply(lambda tmp: convert_to_lat(tmp[1]), axis=1, raw=True)
df['longitude'] = df.apply(lambda tmp: convert_to_long(tmp[0]), axis=1, raw=True)
def handle_data(all_coords, v_wind_rndm, v_plane_rndm, wind_limit=None):
def wind_envelope(df, cutoff):
bool_list = df['v_w'] <= cutoff
return df[bool_list]
# landing coordinates
results_df = pd.DataFrame.from_records([rotate(np.array([m_to_ft(Xi) for Xi in X]),
np.radians(heli_direction) - np.pi/2)
for X in all_coords],
columns=['x', 'y'])
#results_df = pd.DataFrame.from_records([np.array([m_to_ft(Xi) for Xi in X])
# for X in all_coords],
# columns=['x','y'])
polar_coord(results_df)
# plane velocities
results_df['v_plane'] = [ms_to_knots(vel) for vel in v_plane_rndm]
# wind magnitude
results_df['v_w'] = [ms_to_knots(vel[0]) for vel in v_wind_rndm]
# wind direction (is this on unit circle? need to think more...)
results_df['deg_w'] = [vel[1] for vel in v_wind_rndm]
lat_long(results_df)
if wind_limit != None:
print("Cutting out all runs with winds above "+str(wind_limit)+" m/s.")
results_df = wind_envelope(results_df, wind_limit)
print(results_df.describe())
return results_df
def translate_cluster(df, offset):
df['x'] = df['x'].apply(lambda x: x + offset[0])
df['y'] = df['y'].apply(lambda y: y + offset[1])
polar_coord(df)
lat_long(df)
# i don't want to talk about how ugly this code block is lol
def rotate_cluster(df, theta):
array = df.loc[:, ['x', 'y']]
X = [x for x in array['x']]
Y = [y for y in array['y']]
new = [rotate(np.array([x, Y[i]]), theta) for i, x in enumerate(X)]
X = [x[0] for x in new]
Y = [y[1] for y in new]
df['x'] = X
df['y'] = Y
polar_coord(df)
lat_long(df)
def get_mean(df):
mean_x = df.apply(lambda pos: pos[2] * np.cos(pos[3]), axis=1, raw=True).describe()['mean']
mean_y = df.apply(lambda pos: pos[2] * np.sin(pos[3]), axis=1, raw=True).describe()['mean']
return mean_x, mean_y
In [5]:
def visualize(results_df):
# these are not the means of the coordinates, taken individually
# this is the coordinate of the mean impact point, hence the polar coordinates
mean_x, mean_y = get_mean(results_df)
print('Mean Impact: (', mean_x, ',', mean_y, ')')
# assume that the bias is 0 until we have an actual target point
# so CEP = sqrt(MSE) = sqrt(var_norm + var_theta + cov_norm_theta)
# this is the mean distance from the mean impact point
variance = results_df.var(axis=0)[2:4].sum()
covariance = results_df.cov().iloc[2, 3]
# GPS is 95% inside a 2x2 m box, and 95% of values on one axis are within 2*RMS of ideal
# thus RMS = 1/2, since 2* 1/2 =1, and 1 on each side gives a 2 m bound. Therefore MSE = 1/4, which we double
# because there are two axes to account for
GPS_MSE = m_to_ft(m_to_ft(1/4))
bias = 0
total_MSE = variance + covariance + bias + 2*GPS_MSE
CEP = np.sqrt(total_MSE)
print('CEP:', CEP.round(3), 'ft')
plt.figure(1)
plt.hist2d(results_df['v_w'], results_df['deg_w'], 25)
plt.title('wind vectors')
plt.xlabel('v (knots)')
plt.ylabel('direction (degrees)')
plt.show()
plt.figure(2)
plt.hist(results_df['v_w'], 30)
plt.xlabel('wind magnitude (knots)')
plt.ylabel('count')
plt.title('wind speed distribution')
plt.show()
plt.figure(3)
plt.hist(results_df['norm'], 30)
plt.title('distances from origin')
plt.xlabel('radius (ft)')
plt.ylabel('count')
plt.show()
plt.figure(4)
plt.hist2d(results_df['x'],results_df['y'], 30)
plt.plot(mean_x, mean_y, 'r^')
plt.title('landing sites')
plt.xlabel('x (ft)')
plt.ylabel('y (ft)')
plt.show()
plt.figure(5)
results_df.plot.scatter('x', 'y', s=1, c='v_plane', colormap='inferno')
plt.axis('equal')
circle1 = plt.Circle((mean_x, mean_y), CEP, color='r', fill=False)
circle2 = plt.Circle((mean_x, mean_y), 2*CEP, color='tab:orange', fill=False)
circle3 = plt.Circle((mean_x, mean_y), 3*CEP, color='y', fill=False)
plt.gca().add_artist(circle1)
plt.gca().add_artist(circle2)
plt.gca().add_artist(circle3)
plt.title('impact points (and initial helicopter velocities)')
plt.show()
plt.figure(6)
results_df.plot.scatter('x', 'y', s=1, c='v_w', colormap='inferno')
plt.axis('equal')
circle1 = plt.Circle((mean_x, mean_y), CEP, color='r', fill=False)
circle2 = plt.Circle((mean_x, mean_y), 2*CEP, color='tab:orange', fill=False)
circle3 = plt.Circle((mean_x, mean_y), 3*CEP, color='y', fill=False)
plt.gca().add_artist(circle1)
plt.gca().add_artist(circle2)
plt.gca().add_artist(circle3)
plt.title('impact points (and initial wind velocities)')
plt.show()
return ((mean_x, mean_y), CEP)
In [6]:
def plot_on_map(df, CEP):
lat_0 = 45.208778
long_0 = -123.138722
mean_x, mean_y = get_mean(df)
mean_long = convert_to_long(mean_x) + long_0
mean_lat = convert_to_lat(mean_y) + lat_0
avg_CEP_est = (convert_to_lat(CEP) + convert_to_long(CEP)) / 2
map_options = GMapOptions(lat=lat_0, lng=long_0, map_type="satellite", zoom=16)
plot = GMapPlot(x_range=Range1d(), y_range=Range1d(), map_options=map_options)
plot.title.text = "Impact Points"
plot.api_key = "AIzaSyCYsl3277rp6nR3BUUx2q3Z9R4fyI2J1qo"
#this block doesn't work :(
circle1 = Circle(x=mean_long, y=mean_lat, radius=avg_CEP_est, line_width=2, line_color='red')
circle2 = Circle(x=mean_long, y=mean_lat, radius=avg_CEP_est*2, line_width=2, line_color='green')
circle3 = Circle(x=mean_long, y=mean_lat, radius=avg_CEP_est*3, line_width=2, line_color='blue')
plot.add_glyph(circle1)
plot.add_glyph(circle2)
plot.add_glyph(circle3)
source = ColumnDataSource(
data=dict(
lat=df.latitude.apply(lambda y: y + lat_0).tolist(),
lon=df.longitude.apply(lambda x: x + long_0).tolist(),
v_h=df.v_plane.tolist(),
v_w=df.v_w.tolist()))
color_mapper = LinearColorMapper(palette=Viridis5)
circle = Circle(x="lon", y="lat", size="v_w", fill_color={'field': 'v_h', 'transform': color_mapper},
fill_alpha=0.5, line_color=None)
plot.add_glyph(source, circle)
color_bar = ColorBar(color_mapper=color_mapper, ticker=BasicTicker(),
label_standoff=12, border_line_color=None, location=(0,0))
plot.add_layout(color_bar, 'right')
plot.add_tools(PanTool(), WheelZoomTool())
output_notebook()
show(plot)
In [7]:
amount = 1150
# mission success
succ_coords, succ_wind, succ_plane = mc_init(num_sims=amount, dt=0.05, airspeed=15.43334)
# main chute failure
coords_A, wind_A, plane_A = mc_init(num_sims=amount, dt=0.05, mode='A', airspeed=15.43334)
# drogue failure
coords_B, wind_B, plane_B = mc_init(num_sims=amount, dt=0.05, mode='B', airspeed=15.43334)
# early deployment
coords_C, wind_C, plane_C = mc_init(num_sims=amount, dt=0.05, mode='C', airspeed=15.43334)
In [8]:
succ_data = handle_data(succ_coords, succ_wind, succ_plane, wind_limit=None)
succ_pos, succ_CEP = visualize(succ_data)
plot_on_map(succ_data, succ_CEP)
In [9]:
#rotate_cluster(succ_data, np.pi)
#translate_cluster(succ_data, (-500, -500))
#succ_pos, succ_CEP = visualize(succ_data)
#plot_on_map(succ_data, succ_CEP)
In [10]:
data_A = handle_data(coords_A, wind_A, plane_A, wind_limit=None)
pos_A, CEP_A = visualize(data_A)
plot_on_map(data_A, CEP_A)
In [11]:
data_B = handle_data(coords_B, wind_B, plane_B, wind_limit=None)
pos_B, CEP_B = visualize(data_B)
plot_on_map(data_B, CEP_B)
In [12]:
data_C = handle_data(coords_C, wind_C, plane_C, wind_limit=None)
pos_C, CEP_C = visualize(data_C)
plot_on_map(data_C, CEP_C)
In [13]:
succ_data.to_csv(path_or_buf='./drop_test_sample_data/succ_data.csv')
data_A.to_csv(path_or_buf='./drop_test_sample_data/data_A.csv')
data_B.to_csv(path_or_buf='./drop_test_sample_data/data_B.csv')
data_C.to_csv(path_or_buf='./drop_test_sample_data/data_C.csv')
In [14]:
aggregate = pd.concat(map(pd.read_csv, ['./drop_test_sample_data/succ_data.csv',
'./drop_test_sample_data/data_A.csv',
'./drop_test_sample_data/data_B.csv',
'./drop_test_sample_data/data_C.csv']))
In [15]:
aggregate.drop('Unnamed: 0', axis='columns', inplace=True)
print(aggregate.describe())
In [16]:
agg, agg_cep = visualize(aggregate)
In [ ]: