In [ ]:
import numpy as np
import json
import matplotlib.pyplot as plt
Section 1 - Data Loading and Parsing¶
In [ ]:
# load satellite ephemeris data from the satellite.json file
with open('/content/drive/My Drive/Planet/Optim/satellites.json') as f: data_f = json.load(f)
# load sun ephemeris data from the sun.json file
with open('/content/drive/My Drive/Planet/Optim/sun.json') as f: data_sun = json.load(f)
# load city location data from the targets.json file
with open('/content/drive/My Drive/Planet/Optim/targets.json') as f: data_targets = json.load(f)
In [ ]:
# dictionary labels for satellite ephemeris
dict_def_satellite = ["RV_ecef_km","t_ISO","timestamp","name","timestamp_minute"]
# what we want = [time_minute, x, y, z, vx, vy, vz]
# vector of [7,] for 1441 times for 15 satellites
In [ ]:
# store satellite state data into an array
satellite_state = np.zeros([7,len(data_f[0]['timestamp_minute']),len(data_f)])
t_iso_list = data_f[0]['t_ISO']
sat_name_list = []
# len(data_f[0]['timestamp_minute']) = 1441 measurements
for i in range(len(data_f)):
satellite_state[0,:,i] = np.array(data_f[i]['timestamp_minute'])
satellite_state[1:7,:,i] = np.array(data_f[i]['RV_ecef_km']).T
sat_name_list.append(data_f[i]['name'])
In [ ]:
# check format for a single state for a single timestep for 1 satellite
# ----------------------------------------------------------------------------
print('State:',satellite_state[:,0,0])
print('Shape:',satellite_state[:,0,0].shape)
In [ ]:
# dictionary labels for satellite ephemeris
dict_def_sun = ["RV_ecef_km","t_ISO","timestamp","name","timestamp_minute"]
# what we want = [time_minute, x, y, z, vx, vy, vz]
# vector of [7,] for 1441 times
In [ ]:
# store sun ephemeris data into an array
sun_state = np.zeros([7,len(data_sun['timestamp_minute'])])
sun_state[0,:] = np.array(data_sun['timestamp_minute'])
sun_state[1:7,:] = np.array(data_sun['RV_ecef_km']).T
In [ ]:
# check format for a single timestep for sun
# ----------------------------------------------------------------------------
print('State:',sun_state[:,0])
print('Shape:',sun_state[:,0].shape)
In [ ]:
# dictionary labels for targets
dict_def_targets = ["name","country","population","lat_deg","lon_deg","alt_m","R_ecef_km","id"]
# what we want = [id, population, lat_deg, long_deg, alt_m, x, y, z]
# vector of [8,] for 5000 targets
In [ ]:
targets_array = np.zeros([8,len(data_targets)])
target_id_list = []
for i in range(len(data_targets)):
targets_array[0,i] = np.array(data_targets[i]["id"])
targets_array[1,i] = np.array(data_targets[i]["population"])
targets_array[2,i] = np.array(data_targets[i]["lat_deg"])
targets_array[3,i] = np.array(data_targets[i]["lon_deg"])
targets_array[4,i] = np.array(data_targets[i]["alt_m"])
targets_array[5:8,i] = np.array(data_targets[i]["R_ecef_km"])
target_id_list.append(data_targets[i]["id"])
In [ ]:
# check format for a single target
# ----------------------------------------------------------------------------
print('Target:',targets_array[:,0])
print('Shape:',targets_array[:,0].shape)
In [ ]:
# Looks like targets are repeated in the file
counts, bins, bars = plt.hist(targets_array[0,:],bins=5000)
In [ ]:
# see if all locations in the file are unique
id_list = targets_array[0,:].tolist()
from collections import Counter
a = Counter(id_list)
print(a)
# function to delete repeated targets and get a unique target array
def unique(list1):
unique_list = []
# traverse for all elements
for x in list1:
# check if exists in unique_list or not
if x not in unique_list:
unique_list.append(x)
return unique_list
id_list = targets_array[0,:].tolist()
unique_list = unique(id_list)
num_uniq = len(unique_list)
print('number of unique locations:',num_uniq)
unique_targets = np.zeros([8,num_uniq])
for i in range(num_uniq):
index = np.where(targets_array[0,:]==unique_list[i])
ind_val = index[0][0]
unique_targets[0,i] = np.array(data_targets[ind_val]["id"])
unique_targets[1,i] = np.array(data_targets[ind_val]["population"])
unique_targets[2,i] = np.array(data_targets[ind_val]["lat_deg"])
unique_targets[3,i] = np.array(data_targets[ind_val]["lon_deg"])
unique_targets[4,i] = np.array(data_targets[ind_val]["alt_m"])
unique_targets[5:8,i] = np.array(data_targets[ind_val]["R_ecef_km"])
In [ ]:
num_uniq
Out[ ]:
In [ ]:
# satellite array -> [time_minute, x, y, z, vx, vy, vz]' = 7 X 1441 X 15
print("Sat array shape:",satellite_state.shape)
# sun array -> [time_minute, x, y, z, vx, vy, vz]' = 7 X 1441
print("Sun array shape:",sun_state.shape)
# target array -> [id, population, lat_deg, long_deg, alt_m, x, y, z] = 8 X 5000
print("Unique targets array shape:",unique_targets.shape)
Section 2 - Creating Lookup Tables¶
In [ ]:
# create a cmobined 3D vector for Sat1,....Sat15, Sun [16 X 1441] X 3 (x, y, z) coordinates
coordinates_all = np.zeros([16,1441,3])
for j in range(3):
coordinates_all[-1,:,j] = sun_state[j+1,:]
for i in range(15):
coordinates_all[i,:,j] = satellite_state[j+1,:,i]
In [ ]:
import time
start = time.time()
# create an elevation and azimuth table -> 16 X 1441 X 3868
elevation_table = np.zeros([16,1441,num_uniq])
azimuth_table = np.zeros([16,1441,num_uniq])
slantR_table = np.zeros([16,1441,num_uniq])
# compute values for each location and each time step
for i in range(num_uniq):
cos_phi, sin_phi = np.cos(np.deg2rad(unique_targets[2,i])), np.sin(np.deg2rad(unique_targets[2,i]))
cos_lam, sin_lam = np.cos(np.deg2rad(unique_targets[3,i])), np.sin(np.deg2rad(unique_targets[3,i]))
u, v, w = coordinates_all[:,:,0] - unique_targets[-3,i], coordinates_all[:,:,1] - unique_targets[-2,i], coordinates_all[:,:,2] - unique_targets[-1,i]
t = cos_lam*u + sin_lam*v
xEast = -sin_lam*u + cos_lam*v
zUp = cos_phi*t + sin_phi*w
yNorth = -sin_phi*t + cos_phi*w
r = np.sqrt(xEast**2+yNorth**2)
slantR_table[:,:,i] = np.sqrt(r**2+zUp**2)
elevation_table[:,:,i] = np.rad2deg(np.arctan2(zUp,r))
azimuth_table[:,:,i] = np.mod(np.rad2deg(np.arctan2(xEast,yNorth)), 360)
print("time taken to compute:",time.time()-start)
In [ ]:
print("At t = 0, elevation, azimuth and Range from location id=441 for all satellites + sun")
print("Elevation: ",elevation_table[:,0,0], " \n Azimuth: ", azimuth_table[:,0,0], "\n Range: " ,slantR_table[:,0,0])
# each with shape [16 X 1441 X 5000]
In [ ]:
num_sats = 2
plt.figure(1)
for i in range(num_sats): plt.plot(satellite_state[0,:,0],elevation_table[i,:,0],alpha=0.3)
plt.title("Elevation change for Location_Id=441 with time")
plt.xlabel("timestamp_minute")
plt.ylabel("Elevation angle (deg)")
plt.show()
plt.figure(2)
for i in range(num_sats): plt.plot(satellite_state[0,:,0],azimuth_table[i,:,0],alpha=0.3)
plt.title("Azimuth change for Location_Id=441 with time")
plt.xlabel("timestamp_minute")
plt.ylabel("Azimuth angle (deg)")
plt.show()
plt.figure(3)
for i in range(num_sats): plt.plot(satellite_state[0,:,0],slantR_table[i,:,0],alpha=0.3)
plt.title("Range change for Location_Id=441 with time")
plt.xlabel("timestamp_minute")
plt.ylabel("Range angle (km)")
plt.show()
In [ ]:
# check elevation and azimuth tables with results provided
print("timestamp_minute:",satellite_state[0,364,0])
print("satellite_azimuth_deg:", azimuth_table[3,364,0])
print("satellite_elevation_deg:", elevation_table[3,364,0])
print("sun_azimuth_deg:", azimuth_table[-1,364,0])
print("sun_elevation_deg:", elevation_table[-1,364,0])
Section 3 - Data Visualization and Manipulation¶
In [ ]:
# elevation look up table = [15 sats + sun X 1441 timesteps X 3868 locations]
print('Elevation table shape:',elevation_table.shape)
# azimuth look up table = [15 sats + sun X 1441 timesteps X 3868 locations]
print('Azimuth table shape:',azimuth_table.shape)
# range look up table = [15 sats + sun X 1441 timesteps X 3968 locations]
print('Range table shape:',slantR_table.shape)
In [ ]:
# create array for satellite and sun differently
# satellite tables for elevation and range
sats_elevation_table = np.zeros([15,1441,num_uniq])
sats_elevation_table += elevation_table[0:15,:,:] # [15 X 1441 X 3868]
sats_range_table = slantR_table[0:15,:,:]
# sun table for elevation
sun_elevation_table = np.zeros([1,1441,num_uniq])
sun_elevation_table += elevation_table[-1,:,:] # [1 X 1441 X 3868]
# time vactor
time_vector = satellite_state[0,:,0]
In [ ]:
# modify the sun elevation table to set elevation values to 0 for elev < 5
modif_sun_ele_table = sun_elevation_table[:,:]
modif_sun_ele_table[sun_elevation_table<5]=0
print("Shape of new sun elevation table:",modif_sun_ele_table.shape)
In [ ]:
# modify each table to meet minimum requirement for elevation
modif_sat_ele_table_45 = sats_elevation_table[:,:,:]
modif_sat_ele_table_45[sats_elevation_table<45]=0
# sanity check
print("number of value with elevation > 45:",np.count_nonzero(modif_sat_ele_table_45[:,:,:]>45))
print("Shape of new sat elevation table:",modif_sat_ele_table_45.shape)
In [ ]:
# visualize the elevation for three different locations for the entire time history for all satellites
num_sats = 15
for j in range(3):
plt.figure(j)
for i in range(num_sats): plt.plot(time_vector,modif_sat_ele_table_45[i,:,j],':^',alpha=0.3)
plt.plot(time_vector,modif_sun_ele_table[0,:,j],'--k',alpha=0.3,label='sun elevation')
plt.xlabel("timestamp_minute")
plt.ylabel("Elevation angle (deg)")
plt.legend()
plt.show()
As seen above, even though the elevation is greater 45 deg for location 0 for many of the satellites, it will not always recieve a sun elevation greater than 5 degrees, so we can do the following: Set all the sun elevation angles = 1 where sun elevation for a place is greater than 5 deg and set everything else to 0. Multiply this vector element-wise with the all satellites for all locations, so only the elevations with good sunlight will have non-zero values.
In [ ]:
# set non-zero values to 1 for sun elevation
final_sun_ele = np.zeros_like(modif_sun_ele_table)
final_sun_ele[modif_sun_ele_table != 0] = 1
final_sat_ele = np.zeros_like(modif_sat_ele_table_45)
for i in range(num_uniq):
final_sat_ele[:,:,i] = modif_sat_ele_table_45[:,:,i] * np.tile(final_sun_ele[0,:,i],(15,1))
In [ ]:
# visualize the elevation for three different locations for the entire time history for all satellites
num_sats = 15
for j in range(15):
plt.figure(j)
for i in range(num_sats): plt.plot(time_vector,final_sat_ele[i,:,j],':^',alpha=0.3)
plt.plot(time_vector,90*final_sun_ele[0,:,j],'-k',alpha=0.3,label='sun elevation interval')
plt.xlabel("timestamp_minute")
plt.ylabel("Elevation angle (deg)")
plt.legend()
plt.show()
# can actually verify with plotting only sat index = 3, for id=441, plt.plot(time_vector,final_sat_ele[3,:,j],':^',alpha=0.3) matches with example provided
In [ ]:
count_ev = np.zeros(num_uniq)
for i in range(num_uniq):
count_ev[i] = np.count_nonzero(final_sat_ele[:,:,i])
In [ ]:
# print(np.min(count_ev))
# print(np.max(count_ev))
counts, bins, bars = plt.hist(count_ev,bins=50)
In [ ]:
print(counts)
print(np.mean(count_ev))
print(np.max(bins))
From above, we know that there are 123 locations which cannot be imaged whatsoever. Also at an average there are around 11 possible satellites candidates per location.
In [ ]:
# analyzing the imaging choices for each satellite at various timesteps
time_choice = np.zeros([15,1441])
for i in range(15):
for j in range(1441):
time_choice[i, j] = np.count_nonzero(final_sat_ele[i,j,:])
In [ ]:
# number of options to image for all satellite during epoch
plt.figure()
ax = plt.axes(projection='3d')
zline = time_choice
xline = np.arange(15)
yline = time_vector
X, Y = np.meshgrid(xline, yline)
ax.scatter3D(X, Y, zline)
Out[ ]:
In [ ]:
# number of eligible cities to image for all satellites
num_sats = 15
plt.figure(1)
for i in range(num_sats): plt.plot(time_vector, time_choice[i,:],'-',alpha=0.25)
plt.xlabel("time")
plt.ylabel("number of choices")
plt.show()
In [ ]:
a = np.nonzero(final_sat_ele[:,:,0])
print(a[0])
print(a[1])
print('number of valid candidates:',len(a[0]))
for i in range(len(a[0])):
print('sat_ind:',a[0][i],'time_ind:',a[1][i],'elevation',final_sat_ele[a[0][i],a[1][i],0])
a = np.nonzero(final_sat_ele[:,:,1])
print(a[0])
print(a[1])
print('number of valid candidates:',len(a[0]))
for i in range(len(a[0])):
print('sat_ind:',a[0][i],'time_ind:',a[1][i],'elevation',final_sat_ele[a[0][i],a[1][i],1])
Based on this, we will now create a lookup table which has top three elevations for each location.
Section 4 - Optimization¶
In [ ]:
# for every location - chose top num_options(3) elevations and arrange them in a matrix = [15 sats, 1441 timesteps, 3868 loc, (1st, 2nd, 3rd) choices]
def create_rankings(num_options):
multiplier_array = np.zeros([15,1441,num_uniq,num_options])
# create a vector which returns TRUE if locations has candidate sata and FALSE otherwise
options_available = np.zeros(num_uniq)
for i in range(num_uniq):
ind = np.unravel_index(np.argsort(-1*final_sat_ele[:,:,i], axis=None), [15,1441]) # sorting by max elevation
x_index_vec = ind[0][0:num_options]
y_index_vec = ind[1][0:num_options]
counter_options = np.count_nonzero(final_sat_ele[x_index_vec,y_index_vec,i])
if counter_options != 0:
options_available[i] = 1
for j in range(num_options):
# multiplier_array[x_index_vec[j], y_index_vec[j], i, j] = 1 * targets_array[1,i]
multiplier_array[x_index_vec[j], y_index_vec[j], i, j] = 1
return multiplier_array, options_available
In [ ]:
def evaluator(design_vec):
schedule_table = np.zeros([15,1441]) # -> table of location id's
population_table = np.zeros([15,1441]) # -> table of populations
for i in range(num_uniq):
if ((design_vec[i] != 0) & (options_available[i] == 1)):
dummy = multiplier_array[:,:,i,design_vec[i]-1]
index = np.where(dummy!=0) # find index to check if a value is previously chosen
x_ind, y_ind = index[0][0], index[1][0]
if schedule_table[x_ind, y_ind] == 0:
schedule_table[x_ind, y_ind] = unique_targets[0,i]
population_table[x_ind, y_ind] = unique_targets[1,i]
return schedule_table, population_table
def objective(design):
design = design.astype(int)
sch_table, pop_table = evaluator(design)
# return -1*np.count_nonzero(sch_table) # returns number of cities
return -100*np.sum(pop_table)/934184156.0 # returns percentage of total population covered
In [ ]:
test_des = np.array([1]*3868)
# test_des = np.array([1,2]*1934)
# test_des = np.array([1,2,0,3]*967)
import time
num_options = 3
multiplier_array, options_available = create_rankings(num_options)
start = time.time()
time_table, pop_table = evaluator(test_des)
print('precentage of cities pictured:', np.count_nonzero(time_table)*100/3868,'%')
print('percentage population covered: ', np.sum(pop_table)*100/np.sum(unique_targets[1,:]),'%')
print('time taken to evaluate designs:', time.time() - start,'seconds')
In [ ]:
!pip install geneticalgorithm
import numpy as np
from geneticalgorithm import geneticalgorithm as ga
Section 5 - Tests¶
In [ ]:
# HYPER PARAMETERS FOR THE DESIGN SEARCH SPACE
num_options = 2 # has to be atleast 1
multiplier_array, options_available = create_rankings(num_options)
print('Rankings created!')
varbound = np.array([[0,num_options]]*3868)
algorithm_param = {'max_num_iteration': 30, 'population_size': 10,\
'mutation_probability':0.1,'elit_ratio': 0.01,\
'crossover_probability': 0.5, 'parents_portion': 0.3,\
'crossover_type':'uniform', 'max_iteration_without_improv': 5}
model = ga(function=objective, dimension=3868, variable_type='int', variable_boundaries=varbound, algorithm_parameters = algorithm_param)
# run the GA algorithm
start = time.time()
model.run()
end = time.time()
print('Optimum Found!')
schedule, population = evaluator((model.output_dict['variable']).astype(int))
print('Total time taken to evaluate designs:', end - start)
print('Percentage of cities pictured for the optimal design:', np.count_nonzero(schedule)*100/3868,'%')
print('Percentage population covered for the optimal design: ', np.sum(population)*100/np.sum(unique_targets[1,:]),'%')
In [ ]:
# HYPER PARAMETERS FOR THE DESIGN SEARCH SPACE
num_options = 2 # has to be atleast 1
multiplier_array, options_available = create_rankings(num_options)
print('Rankings created!')
varbound = np.array([[0,num_options]]*3868)
algorithm_param = {'max_num_iteration': 25, 'population_size': 40,\
'mutation_probability':0.1,'elit_ratio': 0.01,\
'crossover_probability': 0.5, 'parents_portion': 0.3,\
'crossover_type':'uniform', 'max_iteration_without_improv': 5}
model = ga(function=objective, dimension=3868, variable_type='int', variable_boundaries=varbound, algorithm_parameters = algorithm_param)
# run the GA algorithm
start = time.time()
model.run()
end = time.time()
print('Optimum Found!')
schedule, population = evaluator((model.output_dict['variable']).astype(int))
print('Total time taken to evaluate designs:', end - start)
print('Percentage of cities pictured for the optimal design:', np.count_nonzero(schedule)*100/3868,'%')
print('Percentage population covered for the optimal design: ', np.sum(population)*100/np.sum(unique_targets[1,:]),'%')
In [ ]:
# HYPER PARAMETERS FOR THE DESIGN SEARCH SPACE
num_options = 3 # has to be atleast 1
multiplier_array, options_available = create_rankings(num_options)
print('Rankings created!')
lower_bound = 1
varbound = np.array([[lower_bound,num_options]]*3868)
algorithm_param = {'max_num_iteration': 25, 'population_size': 40,\
'mutation_probability':0.1,'elit_ratio': 0.01,\
'crossover_probability': 0.5, 'parents_portion': 0.3,\
'crossover_type':'uniform', 'max_iteration_without_improv': 5}
model = ga(function=objective, dimension=3868, variable_type='int', variable_boundaries=varbound, algorithm_parameters = algorithm_param)
# run the GA algorithm
start = time.time()
model.run()
end = time.time()
print('Optimum Found!')
schedule, population = evaluator((model.output_dict['variable']).astype(int))
print('Total time taken to evaluate designs:', end - start)
print('Percentage of cities pictured for the optimal design:', np.count_nonzero(schedule)*100/3868,'%')
print('Percentage population covered for the optimal design: ', np.sum(population)*100/np.sum(unique_targets[1,:]),'%')
final_design_run3 = (model.output_dict['variable']).astype(int)
In [ ]:
final_design_run3
Out[ ]:
In [ ]:
# HYPER PARAMETERS FOR THE DESIGN SEARCH SPACE
num_options = 2 # has to be atleast 1
multiplier_array, options_available = create_rankings(num_options)
print('Rankings created!')
lower_bound = 1
varbound = np.array([[lower_bound,num_options]]*3868)
algorithm_param = {'max_num_iteration': 25, 'population_size': 100,\
'mutation_probability':0.1,'elit_ratio': 0.01,\
'crossover_probability': 0.5, 'parents_portion': 0.3,\
'crossover_type':'uniform', 'max_iteration_without_improv': 5}
model = ga(function=objective, dimension=3868, variable_type='int', variable_boundaries=varbound, algorithm_parameters = algorithm_param)
# run the GA algorithm
start = time.time()
model.run()
end = time.time()
print('Optimum Found!')
schedule, population = evaluator((model.output_dict['variable']).astype(int))
print('Total time taken to evaluate designs:', end - start)
print('Percentage of cities pictured for the optimal design:', np.count_nonzero(schedule)*100/3868,'%')
print('Percentage population covered for the optimal design: ', np.sum(population)*100/np.sum(unique_targets[1,:]),'%')
final_design_run4 = (model.output_dict['variable']).astype(int)
In [ ]:
import sys
import numpy
np.set_printoptions(threshold=sys.maxsize)
print(final_design_run4)
In [ ]:
# HYPER PARAMETERS FOR THE DESIGN SEARCH SPACE
num_options = 3 # has to be atleast 1
multiplier_array, options_available = create_rankings(num_options)
print('Rankings created!')
lower_bound = 1
varbound = np.array([[lower_bound,num_options]]*3868)
algorithm_param = {'max_num_iteration': 25, 'population_size': 100,\
'mutation_probability':0.1,'elit_ratio': 0.01,\
'crossover_probability': 0.5, 'parents_portion': 0.3,\
'crossover_type':'uniform', 'max_iteration_without_improv': 5}
model = ga(function=objective, dimension=3868, variable_type='int', variable_boundaries=varbound, algorithm_parameters = algorithm_param)
# run the GA algorithm
start = time.time()
model.run()
end = time.time()
print('Optimum Found!')
schedule, population = evaluator((model.output_dict['variable']).astype(int))
print('Total time taken to evaluate designs:', end - start)
print('Percentage of cities pictured for the optimal design:', np.count_nonzero(schedule)*100/3868,'%')
print('Percentage population covered for the optimal design: ', np.sum(population)*100/np.sum(unique_targets[1,:]),'%')
final_design_run5 = (model.output_dict['variable']).astype(int)
Section 6 - Export Results¶
In [ ]:
import numpy as np
import sys
with open('/content/drive/My Drive/Planet/Optim/results.npy', 'wb') as f:
np.save(f, final_design_run3)
with open('/content/drive/My Drive/Planet/Optim/results.npy', 'rb') as f:
design3 = np.load(f)
np.set_printoptions(threshold=sys.maxsize)
print(design3)
In [ ]:
# export imaging schedule to a .json file
final_design = design3
schedule, population = evaluator(final_design)
export_dict = []
for i in range(num_uniq):
index = np.where(unique_list[i]==schedule)
if len(index[0]) != 0:
sat_num = index[0].astype(int)[0]
time_st = index[1].astype(int)[0]
img_act_duration = 10 + (90-elevation_table[sat_num,time_st,i])
item = {"target_id": unique_targets[0,i], "satellite_name": sat_name_list[sat_num], "target_lat_deg": unique_targets[2,i],\
"target_lon_deg": unique_targets[3,i], "target_alt_m": unique_targets[4,i],\
"t_ISO": t_iso_list[time_st], "timestamp_minute": time_vector[time_st],\
"satellite_azimuth_deg": azimuth_table[sat_num,time_st,i], "satellite_elevation_deg": elevation_table[sat_num,time_st,i],\
"sun_azimuth_deg": azimuth_table[-1,time_st,i], "sun_elevation_deg": elevation_table[-1,time_st,i],\
"activity_duration_s": img_act_duration}
export_dict.append(item)
In [ ]:
with open('/content/drive/My Drive/Planet/Optim/imaging.json', 'w') as fout:
json.dump(export_dict, fout, indent=4)