Probabilistic Programming Module#
Here, we discuss how to use the probabilistic programming module to synthetically generate CDMs.
import kessler
import pyro
import numpy as np
import dsgp4
#we seed everything for reproducibility
pyro.set_rng_seed(10)
#we define the observing instruments
#GNSS first:
gnss_cov_rtn=np.array([1e-9, 1.115849341564346, 0.059309835843067, 1e-9, 1e-9, 1e-9])**2,
instrument_characteristics_gnss={'bias_xyz': np.array([[0., 0., 0.],
[0., 0., 0.]]), 'covariance_rtn': gnss_cov_rtn}
gnss = kessler.GNSS(instrument_characteristics_gnss)
#and then radar:
radar_cov_rtn=np.array([1.9628939405514678, 2.2307686944695706, 0.9660907831563862, 1e-9, 1e-9, 1e-9])**2
instrument_characteristics_radar={'bias_xyz': np.array([[0., 0., 0.],
[0., 0., 0.]]), 'covariance_rtn': radar_cov_rtn}
radar = kessler.Radar(instrument_characteristics_radar)
tles=dsgp4.tle.load('tles_sample_population.txt')
print(f"loaded {len(tles)} TLEs")
loaded 20 TLEs
# tles=dsgp4.tle.load('tles_sample_population.txt')
# tles_filtered=[]
# for tle in tqdm(tles):
# try:
# dsgp4.initialize_tle(tle)
# if tle.apogee_alt()<560*1e3:
# tles_filtered.append(tle)
# except:
# continue
conjunction_model = kessler.model.ConjunctionSimplified(time0=60727.13899462018,
# max_duration_days=7.0,
# time_resolution=600000.0,
# time_upsample_factor=100,
miss_dist_threshold=5000.0,
prior_dict=None,
t_prob_new_obs=0.96,
c_prob_new_obs=0.4,
cdm_update_every_hours=8.0,
mc_samples=100,
mc_upsample_factor=100,
pc_method='MC',
collision_threshold=70,
# likelihood_t_stddev=[371.068006, 0.0999999999, 0.172560879],
# likelihood_c_stddev=[371.068006, 0.0999999999, 0.172560879],
likelihood_time_to_tca_stddev=0.7,
t_observing_instruments=[gnss],
c_observing_instruments=[radar],
tles=tles,)
trace,iterations=conjunction_model.get_conjunction()
/Users/ga00693/Develop/kessler/kessler/model.py:765: UserWarning: To copy construct from a tensor, it is recommended to use sourceTensor.clone().detach() or sourceTensor.clone().detach().requires_grad_(True), rather than torch.tensor(sourceTensor).
pyro.deterministic('conj',torch.tensor(conj))
/Users/ga00693/Develop/kessler/kessler/util.py:85: UserWarning: Using torch.cross without specifying the dim arg is deprecated.
Please either pass the dim explicitly or simply use torch.linalg.cross.
The default value of dim will change to agree with that of linalg.cross in a future release. (Triggered internally at /Users/runner/miniforge3/conda-bld/libtorch_1738206012956/work/aten/src/ATen/native/Cross.cpp:66.)
h_vec = torch.cross(r_vec, v_vec)
After 1170 iterations, generated event with 2 CDMs
trace.nodes['cdms']['infer']['cdms']
[CCSDS_CDM_VERS = 1.0
CREATION_DATE = 2025-02-21T03:20:09.135184
ORIGINATOR = KESSLER_SOFTWARE
MESSAGE_ID = KESSLER_SOFTWARE_ac66059e-20e6-11f0-b2ee-f2a9f71f7e1a
TCA = 2025-02-22T01:17:28.408663
MISS_DISTANCE = 141514.9960370336
RELATIVE_SPEED = 8783.071096663718
RELATIVE_POSITION_R = -3769.8761935982398
RELATIVE_POSITION_T = -110791.09624826182
RELATIVE_POSITION_N = -87963.71484285413
RELATIVE_VELOCITY_R = -144.7709891581174
RELATIVE_VELOCITY_T = -5424.100602146948
RELATIVE_VELOCITY_N = -6906.555719570857
COLLISION_PROBABILITY = 0.0
COLLISION_PROBABILITY_METHOD = MC
OBJECT = OBJECT1
OBJECT_DESIGNATOR = 76126AP
CATALOG_NAME = 9827
OBJECT_NAME = COSMOS 886 DEB
INTERNATIONAL_DESIGNATOR = 76126AP
EPHEMERIS_NAME = NONE
COVARIANCE_METHOD = CALCULATED
MANEUVERABLE = N/A
ORBIT_CENTER = EARTH
REF_FRAME = ITRF
X = 2980.1734541932515
Y = -5181.782698754123
Z = -3422.116689930527
X_DOT = 4.425575967702924
Y_DOT = -1.3655090728792418
Z_DOT = 5.7673902574303915
CR_R = 434.5847746412686
CT_R = 3804.230266317108
CT_T = 41269.012800060416
CN_R = -224.26822649455545
CN_T = -3739.194423008724
CN_N = 4893.388541408624
CRDOT_R = -2.690240978331646
CRDOT_T = -32.45608595843143
CRDOT_N = 3.3748910113664388
CRDOT_RDOT = 0.026616660879604667
CTDOT_R = -0.45587521823581134
CTDOT_T = -3.891974449268055
CTDOT_N = 0.21602245879137436
CTDOT_RDOT = 0.002711690719702106
CTDOT_TDOT = 0.000479432666082355
CNDOT_R = 0.27513941996415164
CNDOT_T = 0.4561951334680935
CNDOT_N = -0.44165624627125066
CNDOT_RDOT = 0.0004400365334853134
CNDOT_TDOT = -0.0003127338855863759
CNDOT_NDOT = 0.004822710541832547
OBJECT = OBJECT2
OBJECT_DESIGNATOR = 65082PT
CATALOG_NAME = 3462
OBJECT_NAME = TITAN 3C TRANSTAGE DEB
INTERNATIONAL_DESIGNATOR = 65082PT
EPHEMERIS_NAME = NONE
COVARIANCE_METHOD = CALCULATED
MANEUVERABLE = N/A
ORBIT_CENTER = EARTH
REF_FRAME = ITRF
X = 2972.4172682119224
Y = -5103.530721008279
Z = -3539.772814553296
X_DOT = 5.826757151768794
Y_DOT = 4.085711393314795
Z_DOT = -0.9752519996834083
CR_R = 15.608895194536064
CT_R = -1831.7454009406745
CT_T = 387066.15903594147
CN_R = 71.5596813096827
CN_T = -6122.2782569260735
CN_N = 5424.209709946298
CRDOT_R = 2.1418685416237584
CRDOT_T = -416.38358827942716
CRDOT_N = 8.17067859275866
CRDOT_RDOT = 0.45361166277724996
CTDOT_R = -0.005166751853750279
CTDOT_T = -0.1424221310036492
CTDOT_N = -0.019218390871709318
CTDOT_RDOT = 1.609E-05
CTDOT_TDOT = 5.202E-06
CNDOT_R = 0.06170571196677281
CNDOT_T = -8.629383890813989
CNDOT_N = 1.1567101080406283
CNDOT_RDOT = 0.010208457878684038
CNDOT_TDOT = -1.756E-05
CNDOT_NDOT = 0.0021928221709182245,
CCSDS_CDM_VERS = 1.0
CREATION_DATE = 2025-02-21T11:20:08.694966
ORIGINATOR = KESSLER_SOFTWARE
MESSAGE_ID = KESSLER_SOFTWARE_adc65ed4-20e6-11f0-b2ee-f2a9f71f7e1a
TCA = 2025-02-22T01:17:28.408663
MISS_DISTANCE = 119296.4071367849
RELATIVE_SPEED = 8781.77597072965
RELATIVE_POSITION_R = -3392.400548990203
RELATIVE_POSITION_T = -79317.283331336
RELATIVE_POSITION_N = -89044.33131422414
RELATIVE_VELOCITY_R = -154.91187193411943
RELATIVE_VELOCITY_T = -5422.481073457308
RELATIVE_VELOCITY_N = -6905.960506693898
COLLISION_PROBABILITY = 0.0
COLLISION_PROBABILITY_METHOD = MC
OBJECT = OBJECT1
OBJECT_DESIGNATOR = 76126AP
CATALOG_NAME = 9827
OBJECT_NAME = COSMOS 886 DEB
INTERNATIONAL_DESIGNATOR = 76126AP
EPHEMERIS_NAME = NONE
COVARIANCE_METHOD = CALCULATED
MANEUVERABLE = N/A
ORBIT_CENTER = EARTH
REF_FRAME = ITRF
X = 2961.416138954747
Y = -5176.42090647565
Z = -3446.5166602514805
X_DOT = 4.441367932633197
Y_DOT = -1.3893739558637068
Z_DOT = 5.749634804852701
CR_R = 418.4789881653135
CT_R = 2833.654900508839
CT_T = 26275.07871526111
CN_R = 235.96450630189815
CN_T = 1379.1837684683846
CN_N = 3824.354403059597
CRDOT_R = -1.629021196870671
CRDOT_T = -18.917640141448356
CRDOT_N = -0.6906003255748716
CRDOT_RDOT = 0.015134294118617029
CTDOT_R = -0.44696031062349284
CTDOT_T = -2.9384291761156875
CTDOT_N = -0.2522263711844156
CTDOT_RDOT = 0.001641663225350465
CTDOT_TDOT = 0.0004784788899478311
CNDOT_R = 0.15032615582479048
CNDOT_T = 0.0059359771572660164
CNDOT_N = -0.3262798823098171
CNDOT_RDOT = 0.0004967254843707441
CNDOT_TDOT = -0.00017283616485068076
CNDOT_NDOT = 0.005746080835028397
OBJECT = OBJECT2
OBJECT_DESIGNATOR = 65082PT
CATALOG_NAME = 3462
OBJECT_NAME = TITAN 3C TRANSTAGE DEB
INTERNATIONAL_DESIGNATOR = 65082PT
EPHEMERIS_NAME = NONE
COVARIANCE_METHOD = CALCULATED
MANEUVERABLE = N/A
ORBIT_CENTER = EARTH
REF_FRAME = ITRF
X = 2973.071864839906
Y = -5103.077176433223
Z = -3539.878719318852
X_DOT = 5.826399482002183
Y_DOT = 4.0863283565475435
Z_DOT = -0.9747953563722136
CR_R = 19.416112754758714
CT_R = -1063.96729664468
CT_T = 122541.70103023999
CN_R = -35.02940219476802
CN_T = 1550.9822755102928
CN_N = 5255.002744568752
CRDOT_R = 1.2925239374152526
CRDOT_T = -126.57004812480902
CRDOT_N = -2.22107658156981
CRDOT_RDOT = 0.1356785119000294
CTDOT_R = -0.010480352857442607
CTDOT_T = 0.30203170115238487
CTDOT_N = 0.03313019924386015
CTDOT_RDOT = -0.00046144760169145884
CTDOT_TDOT = 7.133E-06
CNDOT_R = -0.01550631317382298
CNDOT_T = -0.7484418552935522
CNDOT_N = 1.171583503661754
CNDOT_RDOT = 0.0008530656635604496
CNDOT_TDOT = 1.209E-05
CNDOT_NDOT = 0.0018826631138264438]
#let's save all cdms to kvn files:
#for i in range(len(trace.nodes['cdms']['infer']['cdms'])):
# trace.nodes['cdms']['infer']['cdms'][i].save(f'event2_{i}.kvn')