Estimating a position of a satellite in orbit

Abstract

Estimating the position of a satellite in orbit using a neural network model. The process comprising a training phase and an execution phase. The training phase obtains historical measurement data 501

Description

TECHNICAL FIELD The invention relates to estimating a position of a satellite in orbit In particular, the invention involves training a neural network architecture to approximate equations of motion describing the orbit and to output satellite position as a function of time, where physical laws and externally validated dynamics are used to train the neural network. BACKGROUND It is important to be able to accurately estimate and predict the movement or state of objects in orbit in space. Such objects, i.e. satellites, can include spacecraft undertaking a specific mission. For instance, it is important to be able to estimate a satellite’s position to mitigate risks while in orbit, e.g. controlling the satellite’s path to ensure collisions with debris are avoided. Traditional methods for determining or calculating orbits involve acquiring a series of observations / measurements of an object and applying analytical or statistical methods to estimate position and trajectory. Such statistical methods include batch least squares and Kalman filter techniques. While such methods may be effective for objects in geocentric orbits - particularly in so-called low Earth orbit (LEO) - they face challenges when applied to cislunar space. This is because there is generally more of a scarcity of available observations, as well as significant third-body effects, in cislunar space. It is known to use machine learning methods - in particular, neural network (NN) approaches - to perform orbit determination. A relatively recent approach uses so-called ‘Physics Informed Neural Networks’ (PINN) to integrate deep learning techniques with the fundamental principles of physics, with the aim of modelling and solving complex systems governed by differential equations. This was first described in ‘Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations’, Raissi et al., J. Comp. Phys., vol. 378 (2019), pp 686-707. PINNs leverage the capacity of NNs to approximate functions and discover patterns, while simultaneously enforcing the constraints imposed by physical laws through the architecture of the NN and its loss function. PINNs can predict the behaviour of physical systems even in scenarios where traditional numerical methods may struggle as a result of high computational costs, lack of precise models, and / or incomplete boundary and initial condition data. PINNs allow for known physics to be incorporated directly into the learning process. This has a dual advantage: PINNs can be trained with fewer data points compared with other (traditional) machine learning methods; and PINNs ensure that obtained solutions are physically consistent with the underlying physical laws. This is particularly useful in the field of orbit determination, where equations governing system dynamics are typically well understood, but solving them analytically or numerically for complex geometries and conditions can be difficult. In one previous approach, a PINN was used for orbit determination where the loss function included a term considering a difference between predictions obtained from the NN and from a simulated orbit, and term considering how the NN prediction deviates from the physical laws defined in governing equations for the system. However, there remains a need to provide more accurate methods for determining the position of objects in orbit. It is against this background to which the present invention is set. SUMMARY OF THE INVENTION According to an aspect of the invention there is provided a method of estimating a position of a satellite in orbit. The method comprising performing steps in a training phase. The method comprises obtaining historical measurement data indictive of a measured satellite position over a defined historical time period. The method comprises simulating satellite position over an overall time period that includes the defined historical time period and a defined simulation time period after the historical time period. The simulated satellite position is determined to match the measured satellite position in the historical time period. The method comprises defining a neural network model configured to output predicted satellite position as a function of time. The neural network model comprises at least one hidden layer each having a plurality of activation function nodes. Each activation function node has a respective node weight associated therewith. The method comprises defining a plurality of training points each including an input-output pair of a time point and corresponding measured satellite position from the historical time period. The method comprises training the neural network model using the training points. Training the neural network model comprises adjusting the node weights, and optionally biases, to optimise a defined loss function. The defined loss function includes a data loss term indicative of a difference between predicted satellite position from the neural network model, and a reference satellite position, at each time point of the plurality of training points. The defined loss function includes a physics loss term that comprises: a neural network residual term obtained by evaluating one or more defined equations of motion describing satellite trajectory as a function of time using predicted satellite position output by the neural network model; and a simulation residual term obtained by evaluating the one or more defined equations of motion using the simulated satellite position. The method may comprise, in an execution phase, executing the trained neural network model to estimate the position of the satellite in orbit over time. Optimising the defined loss function may comprise minimising a difference between the neural network residual term and the simulation residual term. The loss function L may be defined as L = Ld + wLp where Ld is the data loss term, Lp is the physics loss term, and w is a loss function weighting term. Optionally, w <1. Further optionally w <0.1. The physics loss term Lp may be given by M 2 i=j where i is the training points, j is a defined set sampled over a defined domain, M is a defined number of points from the defined domain, fGj(t) is the simulation residual term, is the neural network residual term, t is time, and 0 represents the node weights and biases. The reference satellite position may be the measured satellite position. The data loss term Ld may be given by N Ld = “ rM}2 i=l where i is the training points, TV is a total number of the training points, r^t, 0) is a three-dimensional positional vector, including three-dimensional coordinates in space, having values of predicted satellite position obtained from the neural network model, r^t) is the reference satellite position, t is time, and 0 represents the node weights and biases. The position of the satellite may comprise a three-dimensional positional vector including three-dimensional coordinates in space. The one or more defined equations of motion may be a defined system of one or more partial differential equations. The defined system of partial differential equations may describe an Unperturbed Two-Body Problem. The system of partial differential equations may be defined as r / / 1 (¾ + r(t) = 0 where three-dimensional positional vector r(t) is given by / x(0\ r(0 = y(0 W) / where %, y, z represent three-dimensional coordinates in space of the satellite, t is time, du is the second derivative with respect to time, / / is a mass parameter calculated by multiplying the gravitational constant by a mass of a body being orbited by the satellite, and parameter r is given by r = j x2 + y2 + z2. The simulation residual term fG(t) may be given by dtt + rG(0 = / g(O where rc(t) and rG are values of the three-dimensional positional vector r(t) and parameter r, respectively, obtained from the simulated satellite position. The defined system of partial differential equations may describe a Circular Restricted Three-Body Problem. Training the neural network model may comprises adjusting the node weights and biases and evaluating the loss function using the output from the neural network model in an iterative manner until a stop condition is satisfied. Optionally, the stop condition may include one or more of: a defined number of iterations have been performed; a value of the loss function is less than a defined threshold value; and a rate of decrease in a value of the loss function across iterations is less than a defined threshold rate value. The satellite may be in a Keplerian orbit. The satellite may be in a low Earth orbit during the historical time period. The historical time period may cover an initial launch period of the satellite from Earth. The method may comprise measuring, using one or more measurement sensors, data indicative of the position of the satellite during the historical time period to obtain the historical measurement data. The measurement sensors may include one or more of: an optical sensor; a radar sensor; a satellite navigation sensor, optionally GPS; and a radio-telescope. The simulated satellite position may be obtained using a defined satellite trajectory simulation software tool. The simulated satellite position may be obtained using the General Mission Analysis Tool. The method may comprise, in the execution phase, determining a position in space of one or more further objects and, based on the estimated position of the satellite, determining a collision risk for the satellite with one or more of the further objects. The method may comprise, in the execution phase, outputting a control signal to control a trajectory of the satellite based on the estimated position of the satellite. According to another aspect of the invention there is provided a controller for estimating a position of a satellite in orbit. The controller is configured to perform steps in a training phase. The controller is configured to obtain historical measurement data indictive of a measured satellite position over a defined historical time period. The controller is configured to simulate satellite position over an overall time period that includes the defined historical time period and a defined simulation time period after the historical time period. The simulated satellite position is determined to match the measured satellite position in the historical time period. The controller is configured to define a neural network model configured to output predicted satellite position as a function of time. The neural network model comprises at least one hidden layer each having a plurality of activation function nodes. Each activation function node has a respective node weight associated therewith. The controller is configured to define a plurality of training points each including an input-output pair of a time point and corresponding measured satellite position from the historical time period. The controller is configured to train the neural network model using the training points, wherein to train the neural network model the controller is configured to adjust the node weights and biases to optimise a defined loss function. The defined loss function includes a data loss term indicative of a difference between predicted satellite position from the neural network model, and a reference satellite position, at each time point of the plurality of training points. The defined loss function includes a physics loss term that comprises: a neural network residual term obtained by evaluating one or more defined equations of motion describing satellite trajectory as a function of time using predicted satellite position output by the neural network model; and a simulation residual term obtained by evaluating the one or more defined equations of motion using the simulated satellite position. The controller may be configured to, in an execution phase, execute the trained neural network model to estimate the position of the satellite in orbit over time. According to another aspect of the invention there is provided a method of training a neural network model to estimate a position of a satellite in orbit. The method comprises obtaining historical measurement data indictive of a measured satellite position over a defined historical time period. The method comprises simulating satellite position over an overall time period that includes the defined historical time period and a defined simulation time period after the historical time period. The simulated satellite position is determined to match the measured satellite position in the historical time period. The method comprises defining the neural network model to be configured to output predicted satellite position as a function of time. The neural network model comprises at least one hidden layer each having a plurality of activation function nodes. Each activation function node has a respective node weight associated therewith. The method comprises defining a plurality of training points each including an input-output pair of a time point and corresponding measured satellite position from the historical time period. The method comprises training the neural network model using the training points, wherein training the neural network model comprises adjusting the node weights and biases to optimise a defined loss function. The defined loss function includes a data loss term indicative of a difference between predicted satellite position from the neural network model, and a reference satellite position, at each time point of the plurality of training points. The defined loss function includes a physics loss term that comprises: a neural network residual term obtained by evaluating one or more defined equations of motion describing satellite trajectory as a function of time using predicted satellite position output by the neural network model; and a simulation residual term obtained by evaluating the one or more defined equations of motion using the simulated satellite position. According to another aspect of the invention there is provided a computer program product comprising a neural network model trained according to the method defined above. According to another aspect of the invention there is provided a non-transitory, computer readable storage medium storing instructions thereon that, when executed by one or more computer processors, causes the one or more computer processors to perform the method defined above. BRIEF DESCRIPTION OF THE DRAWINGS Examples of the invention will now be described with reference to the accompanying drawings, in which: Figure 1 schematically illustrates a satellite in orbit around the Earth; Figure 2 schematically illustrates a system for training a neural network architecture to be used to estimate a position of the satellite of Figure 1 in orbit; Figure 3(a) shows a plot of estimated position of the satellite of Figure 1 determined according to a previous, known method, and Figure 3(b) shows a plot of estimated position of the satellite of Figure 1 determined according to an example method in accordance with the invention; Figure 4(a) shows a plot of mean absolute error over time of the estimated position determined according to the previous, known method in Figure 3(a) compared against a measured position of the satellite, and Figure 4(b) shows a plot of mean absolute error over time of the estimated position determined according to the example method in Figure 3(b) compared against the measured position of the satellite; and Figure 5 show the steps of a method of determining a position of the satellite of Figure 1 in accordance with the present invention. DETAILED DESCRIPTION Figure 1 schematically illustrates an object 101 in orbit 102 around an origin. In the illustrated example the object 101 is in orbit about the Earth 103. The object I satellite 101 may for instance be a spacecraft. In examples of the invention, the satellite 101 may be regarded as being in cislunar space, i.e. the region of space from the Earth out to and including the region around the surface of the Moon or the Moon’s orbit. In the described example, the satellite 101 is considered to be in Keplerian orbit, i.e. the motion of the satellite 101 relative to the Earth 103 is an ellipse, parabola or hyperbola, forming a two-dimensional orbital plane in three-dimensional space. It is desired to determine I estimate the position of the satellite 101 in space. The state of the satellite 101 may be represented with the position vector r(t), where r(t) = y(t) W) / Here, x,y,z represent the three-dimensional coordinates in space of the satellite 101, and t e (0,T] is time up to an end time T of an orbit under consideration. The dynamics of the satellite 101 may be described by a system of (one or more) differential equations, typically partial differential equations (PDEs). In a first example, the satellite dynamics may be described by the so-called Unperturbed Two-Body Problem (UTBP). This refers to a simplified model in classical mechanics where two bodies interact only through gravitational force, ignoring all other forces or perturbations. This idealized scenario allows for the derivation of exact solutions to their motion, most notably Kepler’s laws of planetary motion. It is applied in situations where the influence of other celestial bodies or forces can be neglected, such as in the initial analysis of satellite orbits. Adopting the Euler’s notation for the partial derivatives the equations of motion for the UTBP can be expressed as the following system of PDEs: K+^-]= 0 where and / / is a mass parameter calculated by multiplying the gravitational constant by the mass of the body that is being orbited, i.e. the Earth 103 in this case. In another example, in order to handle more complex situations the satellite dynamics may be described by the so-called Circular Restricted Three Body Problem (CR3BP). The equations of motion for the CR3BP can be expressed as the following system of PDEs: 1 — jU Az dttx — 2 dty — x -I--5— (x + / / ) -I—rl / 1 — u u d tty + 2 dtx — y + I —— -I--: \ r2; / 1- / / u. \ dttz + z(——+ — = where / / = —.........................., = V(x + A)2 + y2 + z2, r2 = y / (x — 1 + / / )2 + y2 + z2 Here, m1 is the mass of the Earth 103 and m2 is the mass of the Moon (not shown in Figure 1). It will be understood that the satellite dynamics of the object 101 may be described by any suitable system of PDEs. In each case, the system of PDEs needs to be solved for x,y,z to obtain the position of the satellite 101 at a given time. In the present disclosure, a neural network (NN) is trained to solve the appropriate system of PDEs. The NN takes a temporal variable as input and returns the related vector position of the object 101. In this way, the NN acts as a PDE solver while its output represents the approximate solution. The NN can be denoted as r(t, 0), where 0 is a trainable parameter that contains information related to weights and biases. Specifically, in the present disclosure a ‘Physics Informed Neural Network’ (PINN) is trained as a PDE solver. PINNs can extrapolate from (possibly limited) training data while ensuring that the outputs remain consistent with the underlying physics. In particular, this is achieved by incorporating the underlying physics - via the defined equations of motion - into a loss function used to train the PINN. A Physic-Informed Neural Network Approach to Orbit Determination’, Scorsoglio et al., J. Astronautical Sciences, 2023, 70(4), Art. 25 discloses a method that uses a PINN for solving orbit determination problems. A loss function is defined that uses a least-squares estimate as a baseline, and a regularising term is added based on the differential equations modelling the dynamics of the problem. In particular, in Scorsoglio et al. the loss function LSco = MSEU + MSEf has a first mean squared error loss (or data loss) term MSEU and a second mean squared error loss (or physics loss) term MSEf. Specifically, in Scorsoglio et al. the data loss and physics loss terms in the loss function are defined as Nu MSEU = -^^\u(tlu,xlu) - U i = l Nf 1 X-1 2 7 i=i where is the number of training points, Nf is a defined / selected number of points in the entire domain, {t^, xlu, u' denote initial and boundary training data on the solution u(t, x) -* I— 1 r • - x JV f and [tpXf].^ denotes the collocation points within the domain for = ut + N[u] where u is the latent solution and Af[-] is a nonlinear differential operator. During training of the PINN, it is desired adjust the trainable parameters to minimise the loss function. In Scorsoglio et al. the data loss term takes the output of the NN, u^t^x^), and minimises the error relative to a reference ul. In particular, the reference is data obtained by executing a suitable simulation of the orbit of the object under consideration and obtaining the simulated solution at the training points. The present invention is advantageous in that it provides a more accurate method for determining the position of objects in orbit. In particular, the method of the present invention provides a deeper layer of physical accuracy than that provided by a conventional PINN framework (described above). This is achieved by incorporating externally validated dynamics into the PINN, which goes beyond the conventional practice of embedding physical laws solely through differential equations within the NNs architecture. In more detail, the method of the present invention integrates simulated object orbit data into a PINN as an additional term in the loss function. In this way, unlike in a conventional PINN approach, embedding of physical laws in the NN is not reliant only on differential equations within the NN’s architecture. Specifically, the invention beneficially adds terms derived from simulated data into the loss function, and therefore does not rely just on the NN’s output and automatic differentiation for physical laws. These and other advantages of the invention will become apparent in the following description. As mentioned above, examples of the present invention use data obtained from a simulated orbit of the satellite 101 in optimisation of the loss function during training of a PINN that is for estimating satellite position. Typically, a path / trajectory of the satellite 101 is predicted based on a sequence of historical data. The historical data may be real data (true data) that can be extrapolated via simulation. The real data is measurement data indicative of the position of the satellite 101. This measurement data may take different forms. For instance, the measurement data may include range and / or range rate observations using a radio-telescope, e.g. from multiple observatories, which can then be used to determine the state / position of the satellite 101 in a known manner. The measurement data may include angular measurements obtained from known astrometric calculations using optical observations (from one or more optical sensors). Further details on these measurement approaches are discussed in Scorsoglio et al. referenced above. When the satellite 101 is closer to the Earth 103, e.g. geocentric orbits, the measurement data may include Global Positioning System (GPS) data or ground-based observations, e.g. using optical or radar measurements. This may be suitable when the satellite 101 is in low Earth orbit (LEO), i.e. Earth centred orbits with an altitude of less than approximately 2000 kilometres. If the historical / real data to be used to simulate a trajectory of the satellite relates to an initial launch trajectory of the satellite from the Earth, then the historical / real data may be data obtained while the satellite is in LEO, meaning that the measurement techniques suitable for observations in LEO mentioned above may be utilised. The simulation may be performed using any suitable simulation / prediction approach. One software tool that allows for such simulation is the General Mission Analysis Tool (GMAT). GMAT is an open-source software system for space mission design, optimisation and navigation (software.nasa.gov / software / GSC-17177-1). GMAT is used for real-world mission design. The software tool contains models of real-world objects, such as spacecrafts, which are then incorporated into a mission sequence where a user inputs commands to perform mission events. Here, a mission sequence may be regarded as including a path / trajectory of a satellite simulated by the simulation tool being used. GMAT simulates satellite trajectories by solving the equations of motion using numerical integration methods, by considering various perturbative forces, and employing high-precision ephemerides. The numerical integration methods may include a Runge-Kutta method to solve the initial value problem and a Gauss-Jackson integrator to propagate the orbits over long periods. The perturbation models used may include gravitational and non-gravitational perturbations. The force models used may influence a satellite’s motion, such as Central Body Gravity, Third Body Gravity and Relativistic Corrections. The high-precision ephemerides may provide the positions and velocities of planets, moons and other bodies at given times. To propagate the state (position and velocity) of the satellite, the techniques used may include Cowell’s method and / or Encke’s method. Cowell’s method involves direct integration of the equations of motion. Encke’s method is used for perturbation analysis where the motion is split into a reference (unperturbed) trajectory and a perturbation component. In examples of the invention, simulated orbit data of the satellite 101 obtained using GMAT is used to train a PI NN to estimate satellite position. The simulated orbit data is obtained using historical measurement / real positional data of the satellite 101. The mission sequence (positional) data of the satellite 101 in GMAT needs to match the measurement data indicative of satellite position during the (historical) time period that is selected. In one example, during this time period the satellite 101 is in LEO and the measurement data may be obtained by a method suitable for LEO (mentioned above). Exact ephemeris data, e.g. obtained from an ephemeris file including measurements of the satellite 101, is used to construct the measurement / true data (for the historical time period) to be used during the PINN training process. This data may include true positions, velocities and times of the satellite 101 during an historical time period, along with other information about the satellite 101, e.g. mass, drag coefficient, etc. In order to generate simulated positional data for use in the PINN, a mission sequence for a satellite / spacecraft may be set up in GMAT such that an initial position, velocity and time of the satellite matches that of the measurement / true data, and the satellite is constructed to share the same mass, drag coefficient, etc., of the real satellite 101. This ensures that the GMAT initial inputs match those of the true data. The satellite is then propagated along an orbit that GMAT calculates for a desired amount of time, and (simulated / predicted) ephemeris data is produced which includes the (simulated / predicted) positions, velocities and times of the satellite 101 along such an orbit. This ephemeris data will be used during optimisation of the loss function when training the PINN, as described below. As defined above, the NN that is to be used to approximate the relevant system of PDEs describing the motion of the satellite 101 is denoted by r(t, 0), where 0 is a trainable parameter that contains information related to weights and biases. To solve the problem, the NN is trained by minimising a loss function L comprising two terms: one, Ld, that fits the observations (data); and one, Lp, that incorporates the physics by embedding the dynamical equations of the system. The loss function L may be expressed as L = Ld + wLp Here, w is a weighting term that that balances the physics loss term, Lp, in order to guard against a situation in which one term dominates the other, avoiding scale mismatch and ensuring physical consistency. As mentioned above, the equations of motion (e.g. the system of PDEs for the U2BP or CR3BP) governing the satellite 101 are to be embedded into the loss function. When embedding the equations of motion into the loss function of a neural network, this takes the name of a ‘residual’. The term ‘residual’ indeed comes from the notion that in an ideal scenario, the value of this physics-based component of the loss function should approach zero, indicating that the NN’s predictions are in full compliance with the underlying physical laws. This residual essentially measures the discrepancy between the observed behaviour as predicted by the NN and the expected behaviour according to physical principles. In the present invention, two residuals are defined and included in the loss function: one for the NN; and one for the predicted / simulated positional data, e.g. from GMAT. Considering the U2BP, the ‘NN residual’ / (t) is defined as K +    = and the ‘simulated positional data residual’ fG(t) is defined as a 3« +A rG(t) = fG(t) 'G That is, f(t) represents the residual produced by the NN, while fG(t) represents the residual produced with the information extracted from the simulation, e.g. from GMAT, as discussed in greater detail below. The time derivatives in the expression for f (t) may be calculated in any suitable manner. In one example, this is performed using an in-built module in PyTorch. See, for instance, ‘Automatic differentiation in machine learning: A survey’, Baydin et al., February 2018. In this approach, derivatives of numerical functions modelled as graphs can be readily calculated. This is particularly efficient as it bypasses the need for approximation methods, thereby eliminating truncation and round-off errors. Consequently, the PI NN operates as a mesh-free technique, removing the need for a discrete mesh in solving PDEs. In examples of the invention, the two terms of the loss function L may be expressed as N La = _ n(0)2 i=l M 2 i=j Here, i represents the initial training points, j is a set sampled over the entire domain, N is the total number of training points, and M is a selected number of points, e.g. random, from the entire domain. In one illustrative example, N may constitute half of a full orbit of a satellite and M may be a full orbit. That is, the physics loss term Lp includes not only the NN residual but also the simulated positional data residual (or GMAT residual) fG](t). In particular, training the PINN involves minimising a difference between the two residual terms fj(t, 0) and Also, in examples of the invention the data loss term Ld considers a difference between the NN predictions and the real / measurement data. The use of real / measurement data in this way means that it is possible to achieve an accuracy higher than that obtained from a simulation tool such as GMAT (which is different from previous approaches). The loss function of examples of the present invention is therefore different in both its formulation, and the sources of data used, compared to previous approaches. The approach for the LI2BP described above may be similarly applied to the CR3BP, with a focus on the individual components (specified above) rather than the entire vector. This adjustment accounts for the fact that, unlike the two-body scenario, each component of the equation behaves consistently across the system, mirroring the uniformity observed in the two-body problem’s analytical treatment. Figure 2 illustrates an example NN system 20 that may be used in examples of the invention and, in particular, illustrates the process for training the NN (which is a deep NN 201). The feedforward NN architecture for one example is summarised in Table 1 below. SoftPlus activation was used, which is a smooth approximation to the ReLU function, but which is differentiable everywhere. Parameter Value Number of Input Neurons 1 Number of Output Neurons 3 Number of Hidden Layers 10 Neurons per Hidden Layer 128 Activation Function Softplus Optimizer Adam Learning rate 1 x 10'5 w 1 x 10-2 Table 1 To enhance the efficiency of the NN’s optimisation process and minimise the time to convergence, all of the data features were normalised between 0 and 1 using the min-max normalisation 202. This scaling is applied individually to each training example, ensuring uniformity in data representation and facilitating smoother and faster model training. The normalisation may take the form: Xi — min (x) x =------------- max (x) — min (x) The normalised input t' is executed by the NN 201 according to the trainable parameter 0, which is a function of the weights w and biases b, to obtain the output f. The automatic differentiation (AD) 203 is performed, then the data loss and physics loss terms are determined as defined above and used to compose the loss function 204 also defined above. It is then checked whether a stop condition 205 is satisfied. This could involve determining whether the total loss calculated via the loss function is less than a prescribed tolerance value, and whether a prescribed (maximum) number of iterations have been performed. If any of the stop conditions are satisfied then the training process is complete (and the NN can be used in an execution phase). If the stop condition(s) is not satisfied then the calculated total loss is used to update the trainable parameter 0. In particular, the loss function calculates the error between the predicted output and a target. The backpropagation algorithm computes the gradients of the loss function with respect to each weight and bias in the network, namely, the trainable parameter. Using the gradients computed by backpropagation, the optimisation algorithm updates the weights and biases. This process may be repeated in an iterative manner until the stop condition is satisfied. A specific example / experiment is now described in which the predictions of a trained PI NN as described above in accordance with the present invention was compared to actual TLE (two-line element set) data, and high-fidelity data produced by GMAT, for a single orbit of the International Space Station (ISS). As will be apparent from the below, the accuracy of the estimated ISS position according to the PI NN approach of the present invention was greater than that of the GMAT prediction. This demonstrates the effectiveness of the described PINN approach in terms of precision and reliability in orbit determination. It is noted that, in recent studies of PINNs applied to orbit determination, GMAT simulations are taken as a reference for model accuracy. By surpassing the accuracy of GMAT, the present approach of incorporating GMAT-derived data into the PINN system therefore emphasises the improvement in accuracy and performance of the described PINN approach. Figures 3(a) and 3(b) show plots comparing accuracy of estimated orbit position by a neural network when embedded physics is not incorporated (plot 301a in Figure 3(a)) and when embedded physics is incorporated using the described PINN approach (plot 301b in Figure 3(b)), respectively, for the U2BP. Both Figures 3(a) and 3(b) illustrate the training / historical data (obtained from real ISS data) 302 used to train the NNs, the actual orbit 303 (using the real ISS data) beyond the training data, and the GMAT simulation 304 of estimated position beyond the training data. Figure 3(a) illustrates that the estimated position 301a using the NN without embedded physics is unable to capture the intrinsic structure of the system beyond the training domain. On the other hand, Figure 3(b) illustrates that the estimated position 301b when embedded physics is incorporated tracks the dynamics of the system even when exposed to new data that was not part of the training. Moreover, Figure 3(b) illustrates that the GMAT simulation 304 also tracks the actual orbit 303 relatively closely like the PINN estimated position 301b. Figures 4(a) and 4(b) show the mean absolute error along the trajectory of each position component x,y,z of the GMAT simulation and PINN prediction, respectively, against the real ISS data. For both GMAT and PINN, the overall trend is similar. In particular, a sinusoidal pattern is followed over the course of one orbit because of the non-constant velocity of the satellite during its revolution. It is seen that the magnitude of the error is appreciably lower in the case of the PINN estimation in Figure 4(b) compared to GMAT in Figure 4(a). Table 2 below shows a comparison of the overall errors obtained via the PINN and GMAT methods for two different approaches: the LI2BP and the CR3BP. Notably, in both cases the PINN of the present invention demonstrates higher accuracy. Note that, in Scorsoglio et al. mentioned above, GMAT is used as the benchmark. By surpassing GMAT in accuracy, the method of the present invention inherently is also more accurate that the approach of Scorsoglio et al. GMAT error PINN error ISS CR3BP 1.4301194 x IO'3 1.3460317 x IO'3 ISS U2BP 1.4301194 x 10-3 1.1024135 x IO"3 Table 2 Figure 5 shows the steps of a method 50 in accordance with examples of the invention. The method 50 is for estimating a position of a satellite in orbit, e.g. an object in orbit in space about another object, e.g. around the Earth. The method 50 includes computer-implemented steps as described below. At step 501, the method 50 involves obtaining historical measurement data indictive of a measured satellite position over a defined historical time period. The historical time period may be defined to be any suitable amount of time. In one purely illustrative example, the historical time period corresponds to half of one complete orbit of the satellite 101 about the Earth 103. Typically, the historical measurement data is obtained at a plurality of time points in the historical time period. The historical measurement data may be collected using one or more suitable sensors, e.g. optical sensors, radar sensors, radio-telescopes, etc., as mentioned above. The type of sensor used may depend on a position / distance of the satellite 101 relative to the Earth during the historical time period. For instance, in an example in which the satellite 101 is a spacecraft and the historical time period relates to a period immediately following the launch of the spacecraft from Earth, the spacecraft may be in low Earth orbit for the duration of the historical time period such that sensors appropriate for measuring satellite position in LEO, e.g. optical or radar sensors, may be used to obtain the historical measurement data. The historical measurement data may be regarded as true or real satellite position data. At step 502, the method 50 involves simulating the position of the satellite 101 in orbit over an overall time period that includes the defined historical time period and a defined simulation time period. The simulation time period follows the historical time period. The simulation time period may typically be regarded as commencing immediately after the historical time period has finished. The simulated satellite position is obtained using a suitable software tool or program for performing such a simulation, e.g. GMAT. Performing the simulation may involve defining one or more characteristics of the satellite 101 e.g. mass, drag coefficient, etc. The simulated satellite position is determined / fixed to match the measured satellite position in the historical time period. The simulation time period may be any suitable time period. In one purely illustrative example, the simulation time period corresponds to one complete orbit of the satellite 101 about the Earth 103. The simulated satellite position may be determined at a plurality of time points in the overall time period. At step 503, the method 50 involves defining a neural network (NN) model, e.g. a deep NN model, configured to output predicted satellite position as a function of time, i.e. the NN takes time as an input. The satellite position is typically a three-dimensional coordinate position, and so the NN typically has an output layer with three components, one of each coordinate dimension. The input layer may be one-dimensional. The NN has one or more hidden layers each having a plurality of activation function nodes. The activation function may be any suitable activation function known in the art, e.g. a SoftPlus activation as mentioned above. Each activation function node has a respective node weight associated therewith. The NN is to be trained using the observed satellite position. As such, at step 504 the method comprises defining a plurality / set of training points - collectively forming training data - to be used to train the NN. Each training point may correspond to a time point / step at which measured satellite positional data has been obtained / acquired. In such a case, each training point may include an input-output pair of a time point / step and corresponding measured satellite position from the historical time period. The NN is to be trained to estimate satellite position in a manner that ensures accordance with physical laws. In this way, the NN is to be a Physics-Informed Neural Network (PINN) as described above. To achieve this, the NN is to be trained using a loss function that embeds, not only physical laws, but also externally validated dynamics, as described above. Specifically, the loss function is defined to include separate ‘data loss’ and ‘physics loss’ terms. The data loss term is indicative of a difference between predicted satellite position from the NN, and a reference satellite position, at each time point of the plurality of training points. In preferred examples of the invention, the reference satellite position is the measured satellite position, i.e. real / true positional data. In such examples, the orbit may be (considered to be) a Keplerian orbit, which can ensure that measurement / real data is obtainable as required. The data loss term may be defined as the expression for Ld defined above. The physics loss term includes two residual terms. A first of these is a ‘neural network residual’ term. This is obtained by evaluating one or more defined equations of motion describing satellite trajectory as a function of time using predicted satellite position output by the NN. The equations may be a system of PDEs, e.g. such as the ones described above, or any other suitable system of dynamic equations for describing motion of a satellite in orbit. The NN residual term may be determined by evaluating a residual obtained when evaluating the defined equations of motion using the output from the NN at each time step over a domain of interest, e.g. the time covered by the overall time period, such as one complete orbit A second of the residuals included in the physics loss term is a ‘simulation residual’ term. This is obtained by evaluating the one or more defined equations of motion using the simulated satellite position, i.e. using the satellite position obtained at each time step, using a simulation tool such as GMAT, over the domain of interest. The physics loss term may involve a consideration of a difference between the NN residual and simulation residual at each time step. Indeed, the physics loss term may be defined as Lp as defined above, i.e. a summation of the square of the difference between the NN and simulation residuals at each evaluation step. The loss function may be defined as a weighted sum of the data loss and physics loss terms. In some examples, the physics loss carries a lower weighting than the data loss term. For instance, the physics loss term may be weighted at least one, and optionally at least two, orders of magnitude less than the data loss term. Returning to Figure 5, step 505 of the method 50 involves training the neural network model using the training points defined in step 504. Training the NN comprises adjusting the node weights to optimise I minimise the defined loss function in a manner that would be understood by the skilled person. In particular, this step typically involves executing the NN having initial / default node weights (and any other biases) with each of the training points and evaluating the loss function. The gradients of the loss function are then used to adjust the node weights and biases in a known manner, e.g. backpropagation with gradient descent. The NN (with updated weights) is then executed again with the training points and the loss function is evaluated again. This process is repeated in an iterative manner with the aim of optimising I minimising the loss function until a stop condition is reached, e.g. a prescribed number of iterations have been performed, the evaluated loss function value is less than a prescribed threshold value, a rate of decrease in loss function value is below a prescribed threshold rate value. Once the training process / phase has been performed, then at step 506 the method 50 may involve executing the trained NN model to estimate / predict the position of the satellite 101. In particular, the NN takes time as input outputs satellite position, e.g. three-dimensional coordinates. The trained NN may be used to estimate satellite position in real time or to estimate future satellite position at a future time. The estimated satellite position may be used in a number of ways. The estimate may be output to a user. Optionally, at step 507 the method 50 may involve determining a risk of collision of the satellite 101 with one or more further objects, e.g. debris, in space based on the estimated position / trajectory of the satellite 101 and a known / predicted position or trajectory of the further objects. Optionally, the method 50 may involve controlling the motion / trajectory of the satellite 101 based on the estimate satellite position, e.g. if there is a risk of collision determined then to plan / initiate avoidance / correction manoeuvres. Such control may be initiated automatically. Alternatively, an operator may initiate such control, i.e. manual control, based on the estimated satellite position / trajectory and / or collision risk that is output. Controlling the state / trajectory of a satellite / spacecraft can be crucial in the initial stages of a launch, for instance. It can also be important for performing manoeuvres such as docking. 5 Many modifications may be made to the described examples without departing from the scope of the appended claims. The above-described method may be implemented by a controller, which may be in the 10 form of any suitable computing device, for instance one or more functional units or modules implemented on one or more computer processors. Such functional units may be provided by suitable software running on any suitable computing substrate using conventional or customer processors and memory. The one or more functional units may use a common computing substrate (for example, they may run on the same server) or separate 15 substrates, or one or both may themselves be distributed between multiple computing devices. A computer memory may store instructions for performing the methods performed by the controller, and the processor(s) may execute the stored instructions to perform the method.

Claims

1. A method of estimating a position of a satellite in orbit, the method comprising:in a training phase:obtaining historical measurement data indictive of a measured satellite position over a defined historical time period;simulating satellite position over an overall time period that includes the defined historical time period and a defined simulation time period after the historical time period, wherein the simulated satellite position is determined to match the measured satellite position in the historical time period;defining a neural network model configured to output predicted satellite position as a function of time, wherein the neural network model comprises at least one hidden layer each having a plurality of activation function nodes, wherein each activation function node has a respective node weight associated therewith;defining a plurality of training points each including an input-output pair of a time point and corresponding measured satellite position from the historical time period;training the neural network model using the training points, wherein training the neural network model comprises adjusting the node weights to optimise a defined loss function, wherein the defined loss function includes:a data loss term indicative of a difference between predicted satellite position from the neural network model, and a reference satellite position, at each time point of the plurality of training points; anda physics loss term that comprises:a neural network residual term obtained by evaluating one or more defined equations of motion describing satellite trajectory as a function of time using predicted satellite position output by the neural network model; anda simulation residual term obtained by evaluating the one or more defined equations of motion using the simulated satellite position;in an execution phase:executing the trained neural network model to estimate the position of the satellite in orbit over time.

2. A method according to Claim 1, wherein optimising the defined loss function comprises minimising a difference between the neural network residual term and the simulation residual term.

3. A method according to Claim 1 or Claim 2, wherein the loss function L is defined asL = Ld + wLpwhere Ld is the data loss term, Lp is the physics loss term, and w is a loss function weighting term; optionally, wherein w <1; further optionally, wherein w <0.1.

4. A method according to any previous claim, wherein the physics loss term Lp is given by Mi=jwhere i is the training points, j is a defined set sampled over a defined domain, M is a defined number of points from the defined domain, fGj(t) is the simulation residual term, [,(1,0) is the neural network residual term, t is time, and 0 represents the node weights and biases.

5. A method according to any previous claim, wherein the reference satellite position is the measured satellite position.

6. A method according to any previous claim, wherein the data loss term Ld is given by N1 X-1 2Ld ~r^ i=lwhere i is the training points, N is a total number of the training points, r^t, O') is a three-dimensional positional vector, including three-dimensional coordinates in space, having values of predicted satellite position obtained from the neural network model, r^t) is the reference satellite position, t is time, and 0 represents the node weights and biases.

7. A method according to any previous claim, wherein the position of the satellite comprises a three-dimensional positional vector including three-dimensional coordinates in space.

8. A method according to any previous claim, wherein the one or more defined equations of motion is a defined system of one or more partial differential equations.

9. A method according to Claim 8, wherein the defined system of partial differential equations describes an Unperturbed Two-Body Problem.

10. A method according to Claim 9, wherein the system of partial differential equations is defined asK + 77]= 0where three-dimensional positional vector r(t) is given byr(t) = y(t)where x, y, z represent three-dimensional coordinates in space of the satellite, t is time, dtt is the second derivative with respect to time, / z is a mass parameter calculated by multiplying the gravitational constant by a mass of a body being orbited by the satellite, and parameter r is given by11. A method according to Claim 10, wherein the simulation residual term fc(t) is givenby-V’ D 'Gwhere rdt) and rG are values of the three-dimensional positional vector rtf) and parameter r, respectively, obtained from the simulated satellite position.

12. A method according to Claim 8, wherein the defined system of partial differential equations describes a Circular Restricted Three-Body Problem.

13. A method according to any previous claim, wherein training the neural network model comprises adjusting the node weights and biases and evaluating the loss function using the output from the neural network model in an iterative manner until a stop condition is satisfied; optionally, wherein the stop condition includes one or more of: a defined number of iterations have been performed; a value of the loss function is less than a defined threshold value; and a rate of decrease in a value of the loss function across iterations is less than a defined threshold rate value.

14. A method according to any previous claim, wherein the satellite is in a Keplerian orbit.

15. A method according to any previous claim, wherein the satellite is in a low Earth orbit during the historical time period.

16. A method according to any previous claim, wherein the historical time period covers an initial launch period of the satellite from Earth.

17. A method according to any previous claim, the method comprising measuring, using one or more measurement sensors, data indicative of the position of the satellite during the historical time period to obtain the historical measurement data.

18. A method according to Claim 17, wherein the measurement sensors include one or more of: an optical sensor; a radar sensor; a satellite navigation sensor, optionally GPS; and a radio-telescope.

19. A method according to any previous claim, wherein the simulated satellite position is obtained using a defined satellite trajectory simulation software tool.

20. A method according to any previous claim, wherein the simulated satellite position is obtained using the General Mission Analysis Tool.

21. A method according to any previous claim, the method comprising, in the execution phase, determining a position in space of one or more further objects and, based on the estimated position of the satellite, determining a collision risk for the satellite with one or more of the further objects.

22. A method according to any previous claim, the method comprising, in the execution phase, outputting a control signal to control a trajectory of the satellite based on the estimated position of the satellite.

23. A controller for estimating a position of a satellite in orbit, the controller being configured to:in a training phase:obtain historical measurement data indictive of a measured satellite position over a defined historical time period;simulate satellite position over an overall time period that includes the defined historical time period and a defined simulation time period after the historical time period, wherein the simulated satellite position is determined to match the measured satellite position in the historical time period;define a neural network model configured to output predicted satellite position as a function of time, wherein the neural network model comprises at least one hidden layer each having a plurality of activation function nodes, wherein each activation function node has a respective node weight associated therewith;define a plurality of training points each including an input-output pair of a time point and corresponding measured satellite position from the historical time period;train the neural network model using the training points, wherein to train the neural network model the controller is configured to adjust the node weights and biases to optimise a defined loss function, wherein the defined loss function includes:a data loss term indicative of a difference between predicted satellite position from the neural network model, and a reference satellite position, at each time point of the plurality of training points; anda physics loss term that comprises:a neural network residual term obtained by evaluating one or more defined equations of motion describing satellite trajectory as a function of time using predicted satellite position output by the neural network model; anda simulation residual term obtained by evaluating the one or more defined equations of motion using the simulated satellite position;in an execution phase:execute the trained neural network model to estimate the position of the satellite in orbit over time.

24. A method of training a neural network model to estimate a position of a satellite in orbit, the method comprising:obtaining historical measurement data indictive of a measured satellite position over a defined historical time period;simulating satellite position over an overall time period that includes the defined historical time period and a defined simulation time period after the historical time period, wherein the simulated satellite position is determined to match the measured satellite position in the historical time period;defining the neural network model to be configured to output predicted satellite position as a function of time, wherein the neural network model comprises at least one hidden layer each having a plurality of activation function nodes, wherein each activation function node has a respective node weight associated therewith;defining a plurality of training points each including an input-output pair of a time point and corresponding measured satellite position from the historical time period; andtraining the neural network model using the training points, wherein training the neural network model comprises adjusting the node weights and biases to optimise a defined loss function, wherein the defined loss function includes:a data loss term indicative of a difference between predicted satellite position from the neural network model, and a reference satellite position, at each time point of the plurality of training points; anda physics loss term that comprises:a neural network residual term obtained by evaluating one or more defined equations of motion describing satellite trajectory as a function of time using predicted satellite position output by the neural network model; anda simulation residual term obtained by evaluating the one or more defined equations of motion using the simulated satellite position.

25. A computer program product comprising a neural network model trained according to the method of Claim 24.

Images