Best way to calculate probability of collision

charlesc · February 14, 2024, 8:46am

Hi all!

I would like to ask for advice before I try something…
I have two spacecraft orbits for which I would like to calculate the probability of collision. I would like to ask if anyone with expertise in doing this has any tips.

I am considering:

a) generating my own CDMs (as it seems this is what the ProbabilityOfCollision module works with)
b) Using the ProbabilityOfCollision method of the Patera2005 Class for example and simply calculating all the required inputs myself.

Is there a better way of doing this?

For context I have written my own estimation routines so I have access to the covariance matrices and could feasibly propagate them to TCA.

Any help/tips/examples much appreciated!
Thank you!!

Charles

Vincent · February 14, 2024, 9:04am

Hello @charlesc ,

I think generating your own CDMs will be overkill if you already have/can propagate orbits and their associated covariances at/to TCA.

I would suggest creating a probability computing instance (i would recommend Laas2015 or Patera2005) and use this method :

compute(Orbit primaryAtTCA, StateCovariance primaryCovariance, Orbit secondaryAtTCA, StateCovariance secondaryCovariance, double combinedRadius)

or

compute(Orbit primaryAtTCA, StateCovariance primaryCovariance, double primaryRadius, Orbit secondaryAtTCA, StateCovariance secondaryCovariance, double secondaryRadius)

Hope it helps you !

Cheers ,
Vincent

EDIT : In fact you have access to all these methods to compute the probability of collision ShortTermEncounter2DPOCMethod (OREKIT 12.0.1 API)

charlesc · February 14, 2024, 9:09am

Hello @Vincent

Ok perfect. I will try them all then !

Thank you this is very helpful!
Charles

charlesc · February 14, 2024, 11:33am

Hello @Vincent ,

Sorry to bother you again. I have been trying to get a simple example working but I keep getting this error:

----> 2 poc = Patera2005.compute(orbit1, covariance1, radius1, orbit2, covariance2, radius2)
TypeError: descriptor 'compute' for 'Patera2005' objects doesn't apply to a 'CartesianOrbit' object

I have tried with KeplerianOrbit, CartesianOrbit, EquinoctialOrbit and CircularOrbit.

Am I maybe using the wrong kind of Orbit? Should I be using a different kind of class?

My full code snippet is here:


#placeholder states and TCA
state1 = [7000000.0, 0.0, 0.0, 0.0, 7000.0, 0.0]
state2 = [7100000.0, 0.0, 0.0, 0.0, 7000.1, 0.0]
tca = datetime.datetime(2020, 1, 1, 0, 0, 0)

orbit1 = CartesianOrbit(PVCoordinates(Vector3D(float(state1[0]), float(state1[1]), float(state1[2])),
                        Vector3D(float(state1[3]), float(state1[4]), float(state1[5]))),
                        FramesFactory.getEME2000(),
                        datetime_to_absolutedate(tca),
                        Constants.WGS84_EARTH_MU)

orbit2 = CartesianOrbit(PVCoordinates(Vector3D(float(state2[0]), float(state2[1]), float(state2[2])),
                        Vector3D(float(state2[3]), float(state2[4]), float(state2[5]))),
                        FramesFactory.getEME2000(),
                        datetime_to_absolutedate(tca),
                        Constants.WGS84_EARTH_MU)

radius1 = 1.0
radius2 = 1.0

covariance1 = MatrixUtils.createRealMatrix(6, 6)
covariance2 = MatrixUtils.createRealMatrix(6, 6)

#fill covariance with random values
for i in range(6):
    for j in range(6):
        covariance1.setEntry(i, j, random.random())
        covariance2.setEntry(i, j, random.random())

# Patera2005.compute(Orbit primaryAtTCA, StateCovariance primaryCovariance, double primaryRadius, Orbit secondaryAtTCA, StateCovariance secondaryCovariance, double secondaryRadius)
poc = Patera2005.compute(orbit1, covariance1, radius1, orbit2, covariance2, radius2)
print(f"Probability of collision: {poc}")

Thanks in advance!

Charles

EDIT: just realized my covariance is not passed as StateCovariance. Will try fixing that now! Sorry

EDIT2: I get a new error now…

InvalidArgsError: (<class 'org.orekit.ssa.collision.shorttermencounter.probability.twod.Patera2005'>, 'compute', (<CartesianOrbit: Cartesian parameters: {P(7000000.0, 0.0, 0.0), V(0.0, 7000.0, 0.0)}>, <StateCovariance: org.orekit.propagation.StateCovariance@6601cc93>, 1.0, <CartesianOrbit: Cartesian parameters: {P(7100000.0, 0.0, 0.0), V(0.0, 7000.1, 0.0)}>, <StateCovariance: org.orekit.propagation.StateCovariance@63716833>, 1.0))

From this code:

state1 = [7000000.0, 0.0, 0.0, 0.0, 7000.0, 0.0]
state2 = [7100000.0, 0.0, 0.0, 0.0, 7000.1, 0.0]
tca = datetime.datetime(2020, 1, 1, 0, 0, 0)

orbit1 = CartesianOrbit(PVCoordinates(Vector3D(float(state1[0]), float(state1[1]), float(state1[2])),
                        Vector3D(float(state1[3]), float(state1[4]), float(state1[5]))),
                        FramesFactory.getEME2000(),
                        datetime_to_absolutedate(tca),
                        Constants.WGS84_EARTH_MU)

orbit2 = CartesianOrbit(PVCoordinates(Vector3D(float(state2[0]), float(state2[1]), float(state2[2])),
                        Vector3D(float(state2[3]), float(state2[4]), float(state2[5]))),
                        FramesFactory.getEME2000(),
                        datetime_to_absolutedate(tca),
                        Constants.WGS84_EARTH_MU)

radius1 = 1.0
radius2 = 1.0

cov_mat1 = MatrixUtils.createRealMatrix(6, 6)
cov_mat2 = MatrixUtils.createRealMatrix(6, 6)

for i in range(6):

    for j in range(6):
        cov_mat1.setEntry(i, j, random.random())
        cov_mat2.setEntry(i, j, random.random())

covariance1 = StateCovariance(cov_mat1, datetime_to_absolutedate(tca), FramesFactory.getEME2000(), OrbitType.CARTESIAN, PositionAngleType.TRUE)
covariance2 = StateCovariance(cov_mat2, datetime_to_absolutedate(tca), FramesFactory.getEME2000(), OrbitType.CARTESIAN, PositionAngleType.TRUE)

# Patera2005.compute(Orbit primaryAtTCA, StateCovariance primaryCovariance, double primaryRadius, Orbit secondaryAtTCA, StateCovariance secondaryCovariance, double secondaryRadius)
patera2005 = Patera2005() 
poc_result = patera2005.compute(orbit1, covariance1, radius1, orbit2, covariance2, radius2)
print(f"Probability of collision: {poc_result}")

EDIT3:

Fixed it by changing the method call to:

patera2005 = Patera2005() 
poc_result = patera2005.compute(orbit1, covariance1, orbit2, covariance2, radius2, 1e-10)

and putting real values in my covariance matrix (otherwise i was getting some non symmetric matrix error).

Vincent · February 14, 2024, 2:10pm

No worries, is everything ok now ? I apologize for not seeing your post earlier .

charlesc · February 14, 2024, 2:28pm

Yes all working!

Please don’t apologise. You are always so helpful on this forum

Thanks again.
Charles

daiana · February 26, 2024, 7:41am

Hello!! I’m trying to reproduce your code. How did you get the result if ProbabilityOfCollision compute(double xm,
double ym,
double sigmaX,
double sigmaY,
double radius,
org.hipparchus.analysis.integration.UnivariateIntegrator integrator,
int customMaxNbOfEval)
has absolutely different args in comparison with what you used?? Thus I get InvalidArgsError

Vincent · February 28, 2024, 1:37pm

Hi @daiana and sorry for the delay,

If i understand your problem correctly, you are trying to compute a probability of collision using the following signature ?

ProbabilityOfCollision compute(double xm,double ym,double sigmaX,double sigmaY,double radius, UnivariateIntegrator integrator,int customMaxNbOfEval)

This signature is specific to numerical probability of collision computing method and is not necessary strictly speaking (you can also use the default signature which will use default values for the integrator and maximum number of evaluations…)

Could you provide a sample of your code which reproduce the issue ?

Cheers,
Vincent

charlesc · February 28, 2024, 2:20pm

Hi @daiana,

Here is a Gist of the code I used in its entirety: PoC_test.ipynb · GitHub

I’m sorry it is quite messy. It was just a quick Jupyter Notebook to understand how this worked. Hopefully this helps.

Charles

Vincent · March 1, 2024, 8:18am

Hi @charlesc and thank you for sharing this notebook ! I’m sure it will be super helpful to other users .

daiana · March 6, 2024, 11:02am

thank you @Vincent, @charlesc for your response!! All works now. I’m just wondering where did you get your covariance_values3 matrix from?? I’ve tried to get correlation coefficient matrix from BatchLS estimation, but it gives me matrix with diagonal elements 1 and rest are zeroes

charlesc · March 6, 2024, 5:55pm

Hi @daiana,

Glad this was helpful. My covariance values were example values that I had lying around. I don’t have that much experience with the BLS estimation routine in Orekit so I’m not sure I can help you with that. Usually getting an identity matrix suggests “perfect” estimate which does seem suspect… Have you checked out this tutorial by @yzokras. It is a great starting point.

Best,
Charles