Physics-informed Conditional Normalizing Flows for Angles-only Cislunar Orbit Determination

# PHYSICS-INFORMED CONDITIONAL NORMALIZING FLOWS FOR ANGLES-ONLY CISLUNAR ORBIT DETERMINATION
Source: [https://arxiv.org/html/2606.30936](https://arxiv.org/html/2606.30936)
Walther LitteriMassimiliano VasileProfessor and Director, Aerospace Centre of Excellence, University of Strathclyde, G1 1XQ Glasgow, UK

###### Abstract

Generative Astrodynamics is advanced in this work by extending generative modelling to an orbit determination problem in the cislunar environment\. The task is formulated as conditional density estimation, aiming to infer the probability distribution of the initial state from angles\-only measurements over short observation arcs\. A normalising flow is trained on perturbed topocentric observations from Near Rectilinear Halo Orbits, enabling a flexible and potentially multimodal posterior representation\. Given new measurements, the learned density is sampled to generate statistically consistent and physics\-informed state hypotheses\. These estimates are refined via nonlinear least\-squares minimisation, providing a competitive warm start for classical algorithms\.

## 1Introduction

The problem of recovering the orbit followed by a body in space is arguably one of the first scientific problems ever formulated\. With the advent of space exploration and the growing exploitation of the space environment, this problem has also become central to the engineering community, which requires reliable methods to estimate the state of a spacecraft from observational data\. This field of research is known asOrbit Determination\(OD\)\. The literature further distinguishes between different operational scenarios depending on the phase of the mission\[[14](https://arxiv.org/html/2606.30936#bib.bib7)\]\. In particular, when orbit estimation is performed at the very beginning of operations, e\.g\., immediately after launch or following a critical recovery phase, the problem is referred to as Initial Orbit Determination \(IOD\)\. A key difference between standard OD and IOD lies in the availability of prior information\. In routine OD operations, additional data may be available from previously estimated states or measurements provided by onboard instruments\. When such information is unavailable, as in the case of newly launched spacecraft or uncooperative objects, the IOD problem becomes significantly more challenging\.

Over the past decades, interest in space operations in the cislunar environment has increased considerably\. The presence of water ice in permanently shadowed craters near the lunar south pole makes the Moon an attractive location for future scientific and commercial activities\[[1](https://arxiv.org/html/2606.30936#bib.bib8)\]\. The Artemis programme represents the most recent international effort aimed at returning humans to the lunar vicinity, with the crewed lunar flyby mission Artemis II now scheduled for April 2026\[[3](https://arxiv.org/html/2606.30936#bib.bib9)\]\. In the original Artemis architecture, the Gateway lunar orbital outpost was conceived as a strategic counterpart to the Earth‑orbiting International Space Station, designed to support subsequent surface missions, including Artemis IV and V\. However, following a comprehensive restructuring of the Artemis programme announced in early 2026, the role and priority of the Gateway have become significantly less certain, with current plans shifting towards a more direct and sustained human presence on the lunar surface\.

Operating in the cislunar environment introduces several additional challenges compared to Earth\-orbiting missions\. The large Earth–Moon distance complicates communication and synchronisation, introducing significant time delays between command transmission and reception\. Moreover, the intrinsic dynamical instability and chaotic behaviour that characterise many cislunar trajectories further increase operational complexity, reducing the margin for allowable errors and demanding higher levels of precision and accuracy\. Orbit Determination in cislunar space is further complicated by the absence of a dedicated Global Navigation Satellite System \(GNSS\) capable of supporting navigation in the region, as well as by limited communication opportunities that restrict the availability of frequent tracking measurements\.

Traditional approaches to orbit determination typically rely on correction schemes, such as minimising the nonlinear mean squared error in a statistical framework\. The Extended Kalman Filter \(EKF\) integrates knowledge of the dynamics to improve estimation accuracy\. More recently, machine learning and statistical learning techniques have been explored to model unmodelled forces through nonlinear latent force models, which act as filters for uncorrelated Gaussian noise\[[6](https://arxiv.org/html/2606.30936#bib.bib10),[2](https://arxiv.org/html/2606.30936#bib.bib14)\]\. Physics\-Informed Neural Networks \(PINNs\) are neural networks trained with an additional term enforcing the underlying physical laws\[[16](https://arxiv.org/html/2606.30936#bib.bib12),[17](https://arxiv.org/html/2606.30936#bib.bib13)\]\. In the context of OD, PINNs have been used to approximate orbit propagation and measurement models in the cislunar environment, thereby facilitating the solution of the inverse problem\[[13](https://arxiv.org/html/2606.30936#bib.bib11)\]\.

Generative Artificial Intelligence \(AI\) refers to deep learning strategies used to model probability distributions and to sample from them\. Applications of this framework span a variety of domains, from the generation of text, images, and audio\[[10](https://arxiv.org/html/2606.30936#bib.bib15)\]\. In astrodynamics, within the emerging field of Generative Astrodynamics, generative methods based on Variational Autoencoders have been successfully applied to analyse and generate periodic orbits in the Circular\-Restricted Three\-Body Problem \(CR3BP\)\[[8](https://arxiv.org/html/2606.30936#bib.bib17),[9](https://arxiv.org/html/2606.30936#bib.bib16)\], while Normalizing Flow models have been employed for the identification and characterisation of equilibrium solutions in the N\-body problem\[[18](https://arxiv.org/html/2606.30936#bib.bib18)\]\.

This work presents a novel application of Generative Astrodynamics to the solution of an orbit determination problem in the cislunar environment\. Topocentric states corresponding to an observation arc are modelled as a distribution conditioned on angles\-only observations\. A Normalizing Flow model is trained on this distribution, and sampling conditioned on observations provides the estimated states, which are further refined using traditional methods\. It will be shown how the inclusion of an additional term in the loss, based on the propagation of the trajectory, helps regularise the physics loss, yielding more accurate results\.

The paper is structured as follows\. First, an overview of the problem is provided, defining the quantities of interest and the underlying assumptions\. The model definition follows, detailing architectural choices, trade\-off, and the implementation of the physics\-informed loss\. Finally, the performance of the proposed approach is discussed, including a comparison with traditional methods, and concluding remarks are presented\.

## 2Cislunar Orbit Determination

### 2\.1Scenario

The mission analysis for the Gateway platform considers a Near\-Rectilinear Halo Orbit \(NRHO\) with a period of approximately 7 days, located around the second collinear Euler–Lagrange pointL2L\_\{2\}of the Earth–Moon Circular Restricted Three\-Body Problem \(CR3BP\), characterised by the mass parameterμ=0\.01215\\mu=0\.01215\. The Equations of Motion \(EOM\) written below are expressed in nondimensional units and in the synodic \(rotating\) reference frameO​x​y​zOxyz, centred at the barycentre of the Earth\-Moon system, withxx\-axis pointing the Moon from the Earth, itszz\-axis in the angular momentum direction, and itsyy\-axis completing the right\-hand set\[[12](https://arxiv.org/html/2606.30936#bib.bib19)\]:

x¨\\displaystyle\\ddot\{x\}=2​y˙\+Ux,\\displaystyle=2\\dot\{y\}\+U\_\{x\},\(1a\)y¨\\displaystyle\\ddot\{y\}=−2​x˙\+Uy,\\displaystyle=\-2\\dot\{x\}\+U\_\{y\},\(1b\)z¨\\displaystyle\\ddot\{z\}=Uz\.\\displaystyle=U\_\{z\}\.\(1c\)The state vector𝐱=\[x,y,z,x˙,y˙,z˙\]T∈ℝ6\\mathbf\{x\}=\[x,y,z,\\dot\{x\},\\dot\{y\},\\dot\{z\}\]^\{T\}\\in\\mathbb\{R\}^\{6\}collects, in order, the position and velocity components of the massless particle relative to the barycentre of the system in the non\-inertial synodic frame\. The quantitiesUx,Uy,UzU\_\{x\},U\_\{y\},U\_\{z\}denote the components of the gradient of the pseudo\-potential functionU​\(x,y,z\)U\(x,y,z\)with respect to the position coordinates:

U​\(x,y,z\)=12​\(x2\+y2\)\+1−μr1\+μr2\+12​μ​\(1−μ\)\.U\(x,y,z\)=\\frac\{1\}\{2\}\(x^\{2\}\+y^\{2\}\)\+\\frac\{1\-\\mu\}\{r\_\{1\}\}\+\\frac\{\\mu\}\{r\_\{2\}\}\+\\frac\{1\}\{2\}\\mu\(1\-\\mu\)\.\(2\)The scalar distances between the particle and the two primary bodies are respectivelyr1r\_\{1\}andr2r\_\{2\}:

r1=\(x\+μ\)2\+y2\+z2,r2=\(x−\(1−μ\)\)2\+y2\+z2\.r\_\{1\}=\\sqrt\{\(x\+\\mu\)^\{2\}\+y^\{2\}\+z^\{2\}\},\\ r\_\{2\}=\\sqrt\{\(x\-\(1\-\\mu\)\)^\{2\}\+y^\{2\}\+z^\{2\}\}\.Figures[1](https://arxiv.org/html/2606.30936#S2.F1)present a selection of ten trajectories belonging to theL2L\_\{2\}\-NRHO family\. These orbits are obtained by numerically integrating the EOM \([1](https://arxiv.org/html/2606.30936#S2.E1)\) over one orbital period, which varies within approximately±1\\pm 1hour around the nominal 7\-day period\.

![Refer to caption](https://arxiv.org/html/2606.30936v1/Figures/Figure_90_2.png)\(a\)Synodic reference frame\. The red stars denote the Euler\-Lagrange pointsL1,L2L\_\{1\},L\_\{2\}and the blue star the Moon
![Refer to caption](https://arxiv.org/html/2606.30936v1/Figures/Figure_92.png)\(b\)Geocentric inertial reference frame\. The red line denotes the Moon’s orbit over 1 sidereal month

Figure 1:Selection of 10 orbits belonging to theL2L\_\{2\}\-NRHO family, propagated for 1 lunar sidereal periodIn Figure[1\(a\)](https://arxiv.org/html/2606.30936#S2.F1.sf1), the trajectories are depicted in the synodic reference frame while in Figure[1\(b\)](https://arxiv.org/html/2606.30936#S2.F1.sf2)the same orbits are represented in Earth\-centred inertial reference frame \(ECI\), denotedOXYZO\_\{\\mathrm\{XYZ\}\}\. The transformation between a state vector𝐱\\mathbf\{x\}expressed in the synodic frame and the corresponding state𝐗\\mathbf\{X\}in the geocentric inertial frame is given by:

𝐗=\[𝛀​\(t\)𝟎3×3𝟎3×3𝛀​\(t\)\]​\(𝐋𝐱\+𝐱C−𝐱E\)\.\\mathbf\{X\}=\\left\[\\begin\{array\}\[\]\{cc\}\\boldsymbol\{\\Omega\}\(t\)&\\mathbf\{0\}\_\{3\\times 3\}\\\\ \\mathbf\{0\}\_\{3\\times 3\}&\\boldsymbol\{\\Omega\}\(t\)\\end\{array\}\\right\]\\left\(\\mathbf\{Lx\}\+\\mathbf\{x\}\_\{\\mathrm\{C\}\}\-\\mathbf\{x\}\_\{\\mathrm\{E\}\}\\right\)\.\(3\)The matrix𝐋\\mathbf\{L\}performs the dimensionalisation of the synodic state, using the characteristic length \(L​ULU\) and velocity \(L​U/T​ULU/TU\) units\. The vector𝐱C\\mathbf\{x\}\_\{\\mathrm\{C\}\}accounts for the Coriolis contribution associated with the rotating frame, with angular velocity vector𝝎s=\[0,0,ωs\]T,ωs=1/T​U\\boldsymbol\{\\omega\}\_\{s\}=\[0,0,\\omega\_\{s\}\]^\{T\},\\omega\_\{s\}=\{1\}/\{TU\}, while𝐱E\\mathbf\{x\}\_\{\\mathrm\{E\}\}represents the Earth’s state in the synodic frame:

𝐱C=L​U​\[𝟎3𝝎s×\[x,y,z\]T\],𝐱E=L​U​\[\[−μ,0,0\]T𝝎s×\[−μ,0,0\]T\]\.\\mathbf\{x\}\_\{\\mathrm\{C\}\}=LU\\left\[\\begin\{array\}\[\]\{c\}\\mathbf\{0\}\_\{3\}\\\\ \\boldsymbol\{\\omega\}\_\{s\}\\times\[x,y,z\]^\{T\}\\end\{array\}\\right\],\\ \\mathbf\{x\}\_\{\\mathrm\{E\}\}=LU\\left\[\\begin\{array\}\[\]\{c\}\[\-\\mu,0,0\]^\{T\}\\\\ \\boldsymbol\{\\omega\}\_\{s\}\\times\[\-\\mu,0,0\]^\{T\}\\end\{array\}\\right\]\.The time\-dependent rotation matrix𝛀​\(t\)∈ℝ3×3\\boldsymbol\{\\Omega\}\(t\)\\in\\mathbb\{R\}^\{3\\times 3\}maps vectors from the synodic frame to the ECI frame\. It is constructed according to the classical Euler33–11–33rotation sequence, parametrised by the Moon’s orbital elements evaluated at the reference epocht0t\_\{0\}:

𝛀​\(t\)=𝛀3​\(Ω0\)​𝛀1​\(i0\)​𝛀3​\(ω0\+ν0\+ωs​t\)\.\\boldsymbol\{\\Omega\}\(t\)=\\boldsymbol\{\\Omega\}\_\{3\}\(\\Omega\_\{0\}\)\\boldsymbol\{\\Omega\}\_\{1\}\(i\_\{0\}\)\\boldsymbol\{\\Omega\}\_\{3\}\(\\omega\_\{0\}\+\\nu\_\{0\}\+\\omega\_\{s\}t\)\.\(4\)The anglesΩ0,i0,ω0,ν0\\Omega\_\{0\},i\_\{0\},\\omega\_\{0\},\\nu\_\{0\}denote, respectively, the Right Ascension of the Ascending Node, the inclination, the argument of perigee, and the true anomaly of the Moon at the reference epocht0t\_\{0\}, as obtained from ephemeris data\. For ease of reference, the matrices below define the elementary rotations by an angleα\\alphaabout axis11\(denoted𝛀1\\boldsymbol\{\\Omega\}\_\{1\}\) and axis33\(denoted𝛀3\\boldsymbol\{\\Omega\}\_\{3\}\)\. The termsccandssdenote, respectively, the trigonometric functionscos⁡α,sin⁡α\\cos\\alpha,\\sin\\alpha\.

𝛀1​\(α\)=\[1000c−s0s1\],𝛀3​\(α\)=\[c−s0sc0001\]\.\\boldsymbol\{\\Omega\}\_\{1\}\(\\alpha\)=\\left\[\\begin\{array\}\[\]\{ccc\}1&0&0\\\\ 0&c&\-s\\\\ 0&s&1\\end\{array\}\\right\],\\ \\ \\boldsymbol\{\\Omega\}\_\{3\}\(\\alpha\)=\\left\[\\begin\{array\}\[\]\{ccc\}c&\-s&0\\\\ s&c&0\\\\ 0&0&1\\end\{array\}\\right\]\.The times instantst≥0t\\geq 0are defined relative to the reference onet0t\_\{0\}, i\.e\. instantt=0t=0corresponds to the the actual \(physical\) timet0t\_\{0\}\.

Table 1:Reference quantities for the Orbit Determination ProblemFinally, Table[1](https://arxiv.org/html/2606.30936#S2.T1)reports the reference units defined in section and used in the following of the work\.

### 2\.2Angles\-Only Observables

Consider the Orbit Determination problem depicted in Figure[2](https://arxiv.org/html/2606.30936#S2.F2)\. An observer is at a position on the Earth’s surface, given by vector𝐑o\\mathbf\{R\}\_\{o\}with respect to the geocentric frameOX​Y​ZO\_\{XYZ\}\. An observation arc is given by the set of observables\{𝐲i\}i=1N\\left\\\{\\mathbf\{y\}\_\{i\}\\right\\\}\_\{i=1\}^\{N\}, withNNthe number of observations conducted at successive time instantstit\_\{i\},i=1,2,…,Ni=1,2,\\ldots,N\. Those observables are derived from the topocentric position of the observer\. The goal of an Orbit Determination problem is to recover the Geocentric state𝐗1∈ℝ6\\mathbf\{X\}\_\{1\}\\in\\mathbb\{R\}^\{6\}of the target object at the beginning of the observation arc, i\.e\., at instantt=t1t=t\_\{1\}\. Let be the geocentric states𝐗i\\mathbf\{X\}\_\{i\}belonging to the observation arc and their respective position and velocity components, namely,𝐑i\\mathbf\{R\}\_\{i\}and𝐕i\\mathbf\{V\}\_\{i\}:

𝐗i=\[𝐑i𝐕i\],i=1,2,…,N\.\\mathbf\{X\}\_\{i\}=\\left\[\\begin\{array\}\[\]\{c\}\\mathbf\{R\}\_\{i\}\\\\ \\mathbf\{V\}\_\{i\}\\end\{array\}\\right\],\\ i=1,2,\\ldots,N\.The line\-of\-sight vector𝝆i\\boldsymbol\{\\rho\}\_\{i\}, is defined as the difference between the target’s position vector and the observer’s position, in the same geocentric reference frames:

𝝆i=𝐑i−𝐑o\.\\boldsymbol\{\\rho\}\_\{i\}=\\mathbf\{R\}\_\{i\}\-\\mathbf\{R\}\_\{o\}\.\(5\)The observables considered in this work are angles\-only measurements, namely the azimuth and elevation angles, denoted byθi\\theta\_\{i\}andφi\\varphi\_\{i\}, respectively:

𝐲i=\[θiφi\],i=1​…,N\.\\mathbf\{y\}\_\{i\}=\\left\[\\begin\{array\}\[\]\{c\}\\theta\_\{i\}\\\\ \\varphi\_\{i\}\\end\{array\}\\right\],\\ i=1\\ldots,N\.These quantities describe the direction of𝝆i\\boldsymbol\{\\rho\}\_\{i\}in a local topocentric frame centred at the observer\. Let\(Xi,Yi,Zi\)\(X\_\{i\},Y\_\{i\},Z\_\{i\}\)denote the components of𝝆i\\boldsymbol\{\\rho\}\_\{i\}expressed in this local frame\. The measurement equations at epochtit\_\{i\}are then given by:

θi\\displaystyle\\theta\_\{i\}=arctan⁡\(YiXi\)\+εθ,\\displaystyle=\\arctan\\left\(\\frac\{Y\_\{i\}\}\{X\_\{i\}\}\\right\)\+\\varepsilon\_\{\\theta\},\(6\)φi\\displaystyle\\varphi\_\{i\}=arcsin⁡\(Zi‖𝝆i‖2\)\+εφ\.\\displaystyle=\\arcsin\{\\left\(\\frac\{Z\_\{i\}\}\{\|\|\\boldsymbol\{\\rho\}\_\{i\}\|\|\_\{2\}\}\\right\)\+\\varepsilon\_\{\\varphi\}\}\.\(7\)
![Refer to caption](https://arxiv.org/html/2606.30936v1/Figures/od_framework.png)Figure 2:Orbit Determination problem by Angles\-Only measurements in ECI reference frameThe quantitiesεθ\\varepsilon\_\{\\theta\}andεφ\\varepsilon\_\{\\varphi\}denote the sensor\-dependent measurement errors, here modeled as Gaussian noise with standard deviationsσθ\\sigma\_\{\\theta\}andσφ\\sigma\_\{\\varphi\}, respectively\. Their values are reported in Table[2](https://arxiv.org/html/2606.30936#S2.T2), together with the observer’s information \(namely, latitude, longitude, and altitude\)\. The observer is the Cebreros Ground Station \(near Madrid, Spain\), operated by the European Space Agency \(ESA\)\.

Table 2:Information on the observer’s location and features
### 2\.3Nonlinear Least\-Squares State Estimation

A standard approach to state estimation in the cislunar environment is Nonlinear Least\-Squares Estimation \(NLSE\) based on a Gauss–Newton iterative scheme\. The method minimises the discrepancy between measured and predicted observables over an observation arc\. In the following, only a concise description is provided; the reader is referred to the relevant literature for a comprehensive treatment\. Let

𝝈i=\[arctan⁡\(sin⁡\(θi−θ^i\)cos⁡\(θi−θ^i\)\)φi−φ^i\]\\boldsymbol\{\\sigma\}\_\{i\}=\\left\[\\begin\{array\}\[\]\{c\}\\arctan\\left\(\\dfrac\{\\sin\(\\theta\_\{i\}\-\\hat\{\\theta\}\_\{i\}\)\}\{\\cos\(\\theta\_\{i\}\-\\hat\{\\theta\}\_\{i\}\)\}\\right\)\\\\ \\varphi\_\{i\}\-\\hat\{\\varphi\}\_\{i\}\\end\{array\}\\right\]\(8\)be the residual vector at observation timetit\_\{i\}, defined as the difference between the true observables\(θi,φi\)\(\\theta\_\{i\},\\varphi\_\{i\}\)and the estimated observables\(θ^i,φ^i\)\(\\hat\{\\theta\}\_\{i\},\\hat\{\\varphi\}\_\{i\}\)\. The azimuth residual is expressed using thearctan⁡\(sin/cos\)\\arctan\(\\sin/\\cos\)formulation in order to properly account for angular periodicity and avoid discontinuities at±π\\pm\\pi\. The estimated observables are obtained by applying the noiseless measurement model of Equations \([6](https://arxiv.org/html/2606.30936#S2.E6)\)\-\([7](https://arxiv.org/html/2606.30936#S2.E7)\) from the observer position at the instanttit\_\{i\}to vector𝝆i\{\\boldsymbol\{\\rho\}\}\_\{i\}, this last computed with Equation \([5](https://arxiv.org/html/2606.30936#S2.E5)\) after numerically integrating the EOM \([1](https://arxiv.org/html/2606.30936#S2.E1)\) from the estimated initial state𝐗^1\\hat\{\\mathbf\{X\}\}\_\{1\}, expressed in synodic coordinates by appropriately inverting the transformation given in Equation \([3](https://arxiv.org/html/2606.30936#S2.E3)\)\. The Gauss–Newton update for the initial state estimate at iterationjjis derived in the following, with𝝈∈ℝ2​N\\boldsymbol\{\\sigma\}\\in\\mathbb\{R\}^\{2N\}denoting global residual vector, obtained by stacking all the residuals at any given instant:

𝐗^1j=𝐗^1j−1−λ​𝐉†​𝝈j−1,j=2,…,Nmax\.\\hat\{\\mathbf\{X\}\}\_\{1\}^\{j\}=\\hat\{\\mathbf\{X\}\}\_\{1\}^\{j\-1\}\-\\lambda\\mathbf\{J\}^\{\\dagger\}\\boldsymbol\{\\sigma\}^\{j\-1\},\\qquad j=2,\\ldots,N\_\{\\max\}\.\(9\)Matrix𝐉\\mathbf\{J\}denotes the Jacobian matrix of the global residual vector to the adimensional initial state

𝐉=∂𝝈j∂𝐱^1j−1∈ℝ6×2​N,\\mathbf\{J\}=\\frac\{\\partial\\boldsymbol\{\\sigma\}^\{j\}\}\{\\partial\\hat\{\\mathbf\{x\}\}\_\{1\}^\{j\-1\}\}\\in\\mathbb\{R\}^\{6\\times 2N\},and\(⋅\)†\(\\cdot\)^\{\\dagger\}the inverse matrix;λ=1\\lambda=1is a damping coefficient introduced to improve convergence robustness\. Note that for computational convenience the adimensional version of the initial state is used\. The inverse of the Jacobian is evaluated using a truncated Singular Value Decomposition \(SVD\)\. If

𝐉=𝐔𝐒𝐕\\mathbf\{J\}=\\mathbf\{U\}\\mathbf\{S\}\\mathbf\{V\}\(10\)is the SVD of the Jacobian, only singular values greater than a prescribed threshold \(here1×10−41\\text\{\\times\}\{10\}^\{\-4\}\) are retained\. The pseudoinverse is then constructed as

𝐉†=𝐕T​𝐒∗​𝐔T,\\mathbf\{J\}^\{\\dagger\}=\\mathbf\{V\}^\{\\mathrm\{T\}\}\\mathbf\{S\}^\{\*\}\\mathbf\{U\}^\{\\mathrm\{T\}\},\(11\)where𝐒∗\\mathbf\{S\}^\{\*\}contains the reciprocals of the retained singular values and zeros elsewhere;\(⋅\)T\(\\cdot\)^\{\\mathrm\{T\}\}denotes the transpose matrix\. The iterative process described in Equation \([9](https://arxiv.org/html/2606.30936#S2.E9)\) is repeated until either a maximum number of iterationsNmax=100N\_\{\\mathrm\{max\}\}=100is reached or convergence is achieved\. Convergence is assessed by monitoring the norm of the difference between successive estimates of the initial state, according to

‖𝐱^1j−𝐱^1j−1‖≤1×10−9\.\\left\\\|\\hat\{\\mathbf\{x\}\}\_\{1\}^\{\\,j\}\-\\hat\{\\mathbf\{x\}\}\_\{1\}^\{\\,j\-1\}\\right\\\|\\leq$1\\text\{\\times\}\{10\}^\{\-9\}$\.

### 2\.4Dataset

A dataset of synthetic measurement arcs is generated starting from the target orbits expressed in the ECI frame, as shown in Figure[1\(b\)](https://arxiv.org/html/2606.30936#S2.F1.sf2)\. More specifically, for each of the 10 selected orbits, 100 observation arcs are constructed, each consisting of 10 observations computed using Equation \([5](https://arxiv.org/html/2606.30936#S2.E5)\)\. For every arc, both the initial observation time and the time interval between consecutive observations are randomly selected in order to produce a diverse set of observation conditions\. Observation cadences include measurements taken every 5 minutes \(corresponding to the integration time step\), every 10 minutes, every 15 minutes, and sequences alternating between 5 and 10 minutes\. As a result, the total duration of the observation arcs ranges from a minimum of 50 min to a maximum of 150 min\. To illustrate the observation scenario, Figure[3](https://arxiv.org/html/2606.30936#S2.F3)shows the target orbit \(in black\) together with the Moon’s orbit in the ECI frame\. A representative subset of 10 observation arcs is also displayed as star markers, with colours distinguishing the different arcs\. This visualization highlights the intrinsic difficulty of the orbit determination problem, which is characterised by relatively short observation intervals that challenge accurate estimation of the orbital state\.

![Refer to caption](https://arxiv.org/html/2606.30936v1/Figures/Figure_94_1.png)Figure 3:ECI frame: example a the target orbit \(black\) with a selection of 10 observation arcs\. In red: Moon’s orbitThe dataset of synthetic observations is then constructed by stacking the time instants together with the angles\-only observables and applying the measurement model of Equations \([6](https://arxiv.org/html/2606.30936#S2.E6)\)–\([7](https://arxiv.org/html/2606.30936#S2.E7)\)\. For each observation arc, 10 measurement realisations are generated by sampling from the noise distribution\. We denote by𝐘∈ℝd×3​N\\mathbf\{Y\}\\in\\mathbb\{R\}^\{d\\times 3N\}the resulting dataset of topocentric observations, while the corresponding initial states are stored in𝐗¯1∈ℝd×6\\bar\{\\mathbf\{X\}\}\_\{1\}\\in\\mathbb\{R\}^\{d\\times 6\}\. Thus, for each row indexk≤dk\\leq d:

𝐘k=\[…,ti,θi,φi,…\],𝐗¯1,k=𝐗1,i=1,…,N,k=1,…,d\.\\mathbf\{Y\}\_\{k\}=\[\\ldots,t\_\{i\},\\theta\_\{i\},\\varphi\_\{i\},\\ldots\],\\ \\bar\{\\mathbf\{X\}\}\_\{1,k\}=\\mathbf\{X\}\_\{1\},\\ i=1,\\ldots,N,\\ \\ k=1,\\ldots,d\.Given the procedure described above, the original dataset consists of 10 000 measurement arcs\. However, to represent a more realistic observation scenario, only feasible measurements are retained\. Specifically, an observation arc is considered valid only if the target elevation is positive at every observation within the arc:

φi≥0,i=1,…,N\.\\varphi\_\{i\}\\geq 0,i=1,\\ldots,N\.The resulting dataset has sized=4119d=4119and is subsequently shuffled and split into training \(ntrain=n\_\{\\mathrm\{train\}\}=3119 samples\), validation \(nval=n\_\{\\mathrm\{val\}\}=500 samples\), and test \(ntest=n\_\{\\mathrm\{test\}\}=500 samples\) subsets\.

## 3Generative Astrodynamics

In this Section the generative modelling of the Orbit Determination problem is presented, introducing first the architecture used and then detailing the implementation of the physics\-informed approach\.

### 3\.1Conditional Normalizing Flow Models

Although the Orbit Determination problem is a well\-known inverse problem, it is also meaningful to interpret it as the task of learning a probability distribution from a finite set of samples, that is, as agenerativemodelling problem\[[7](https://arxiv.org/html/2606.30936#bib.bib3)\]\. This perspective enables a variety of advantages and applications, most notably uncertainty quantification of the estimated target state\. More in detail, the target orbit initial state𝐗¯1,k\\bar\{\\mathbf\{X\}\}\_\{1,k\}can be intended to be the realisation of an unknown probability distribution conditioned to the set of observationsp​\(𝐗¯1,k\|𝐘k\)p\(\\bar\{\\mathbf\{X\}\}\_\{1,k\}\|\\mathbf\{Y\}\_\{k\}\), i\.e\.,𝐗¯1,k∼p​\(𝐗¯1,k\|𝐘k\)\\bar\{\\mathbf\{X\}\}\_\{1,k\}\\sim p\(\\bar\{\\mathbf\{X\}\}\_\{1,k\}\|\\mathbf\{Y\}\_\{k\}\)\. Normalising flow models provide an exact and tractable formulation of this conditional data likelihood through a sequence of invertible transformations\[[4](https://arxiv.org/html/2606.30936#bib.bib2),[11](https://arxiv.org/html/2606.30936#bib.bib4)\]\. Let be𝐰∈ℝ6\\mathbf\{w\}\\in\\mathbb\{R\}^\{6\}a random variable with a known and tractable probability density functionp​\(𝐰\)p\(\\mathbf\{w\}\), and letf:ℝ6→ℝ6f:\\mathbb\{R\}^\{6\}\\rightarrow\\mathbb\{R\}^\{6\}the differentiable and invertible mapping conditioned on the observations:

𝐰=f​\(𝐗¯1,k\|𝐘k\),\\mathbf\{w\}=f\(\\bar\{\\mathbf\{X\}\}\_\{1,k\}\|\\mathbf\{Y\}\_\{k\}\),\(12\)so that

𝐗¯1,k=f−1​\(𝐰\|𝐘k\)=g​\(𝐰\|𝐘k\)\\bar\{\\mathbf\{X\}\}\_\{1,k\}=f^\{\-1\}\(\\mathbf\{w\}\|\\mathbf\{Y\}\_\{k\}\)=g\(\\mathbf\{w\}\|\\mathbf\{Y\}\_\{k\}\)\(13\)is the inverse transformation\. By the change\-of\-variables formula, the log\-likelihood of the target distribution can be written as:

logp\(𝐗¯1,k\|𝐘k\)=logp\(𝐰\)\+log\|𝐉\(f\(𝐗¯1,k\|𝐘k\)\)\|,\\log p\(\\bar\{\\mathbf\{X\}\}\_\{1,k\}\|\\mathbf\{Y\}\_\{k\}\)=\\log p\(\\mathbf\{w\}\)\+\\log\\left\|\\mathbf\{J\}\(f\(\\bar\{\\mathbf\{X\}\}\_\{1,k\}\|\\mathbf\{Y\}\_\{k\}\)\)\\right\|,\(14\)with\|𝐉\(f\(𝐗¯1,k\|𝐘k\)\)\|\\left\|\\mathbf\{J\}\(f\(\\bar\{\\mathbf\{X\}\}\_\{1,k\}\|\\mathbf\{Y\}\_\{k\}\)\)\\right\|denoting the determinant of the Jacobian matrix𝐉​\(f\)∈ℝ6×6\\mathbf\{J\}\(f\)\\in\\mathbb\{R\}^\{6\\times 6\}of the mapping evaluated at state vectors𝐗¯1,k\\bar\{\\mathbf\{X\}\}\_\{1,k\}\.

This work adopts a Real Non\-Volume Preserving \(Real\-NVP\) normalising flow\[[5](https://arxiv.org/html/2606.30936#bib.bib5)\], consisting of a composition ofna=16n\_\{a\}=16affine coupling transformations \(layers\)fj,j=1,…,naf\_\{j\},j=1,\\ldots,n\_\{a\},

f=fna∘fna−1∘…∘f2∘f1,f=f\_\{n\_\{a\}\}\\circ f\_\{n\_\{a\}\-1\}\\circ\\ldots\\circ f\_\{2\}\\circ f\_\{1\},each defined as𝐲j=fj​\(𝐱j\|𝐘k\)\\mathbf\{y\}\_\{j\}=f\_\{j\}\(\\mathbf\{x\}\_\{j\}\|\\mathbf\{Y\}\_\{k\}\):

𝐲j=𝐛j⊙𝐱j\+\(1−𝐛j\)⊙\(𝐱j⊙esj​\(𝐛j⊙𝐱j,𝐘k\)\+tj​\(𝐛j⊙𝐱j,𝐘k\)\),\\mathbf\{y\}\_\{j\}=\\mathbf\{b\}\_\{j\}\\odot\\mathbf\{x\}\_\{j\}\+\(1\-\\mathbf\{b\}\_\{j\}\)\\odot\\\!\\left\(\\mathbf\{x\}\_\{j\}\\odot e^\{\\,s\_\{j\}\(\\mathbf\{b\}\_\{j\}\\odot\\mathbf\{x\}\_\{j\},\\mathbf\{Y\}\_\{k\}\)\}\+t\_\{j\}\(\\mathbf\{b\}\_\{j\}\\odot\\mathbf\{x\}\_\{j\},\\mathbf\{Y\}\_\{k\}\)\\right\),\(15\)where𝐛j∈ℝ6\\mathbf\{b\}\_\{j\}\\in\\mathbb\{R\}^\{6\}is a binary masking vector alternatively masking the position, and velocity components and⊙\\odotdenotes the Hadamard \(element\-wise\) product\. At each coupling layer, the masked components of the input,𝐛j⊙𝐱j\\mathbf\{b\}\_\{j\}\\odot\\mathbf\{x\}\_\{j\}, remain unchanged and serve as conditioning variables\. These components are concatenated with the observation vector𝐘k\\mathbf\{Y\}\_\{k\}and jointly provided as input to the neural networks defining the scale functionsj​\(⋅,𝐘k\)s\_\{j\}\(\\cdot,\\mathbf\{Y\}\_\{k\}\)and translation functiontj​\(⋅,𝐘k\)t\_\{j\}\(\\cdot,\\mathbf\{Y\}\_\{k\}\)\. The resulting scale and translation outputs are then applied exclusively to the complementary \(unmasked\) components\(1−𝐛j\)⊙𝐱j\(1\-\\mathbf\{b\}\_\{j\}\)\\odot\\mathbf\{x\}\_\{j\}, while the masked components are passed through the layer without modification\. This structure ensures a tractable form of the Jacobian, whose log\-determinant is simply written as

log⁡\|𝐉​\(fj\)\|=∑i=16\(1−bj,i\)​sj​\(𝐛j⊙𝐱j,𝐘k\)i\.\\log\|\\mathbf\{J\}\(f\_\{j\}\)\|=\\sum\_\{i=1\}^\{6\}\(1\-b\_\{j,i\}\)\\,s\_\{j\}\(\\mathbf\{b\}\_\{j\}\\odot\\mathbf\{x\}\_\{j\},\\mathbf\{Y\}\_\{k\}\)\_\{i\}\.The learnable transformation functions in this work are parametrised by feedforward neural networks composed of three fully connected layers with hidden dimensiondh=256d\_\{h\}=256and LeakyReLU activations \(negative slope 0\.1\)\. To ensure numerical stability and prevent excessive Jacobian determinants during training, the outputs fromsjs\_\{j\}andtjt\_\{j\}on the unmasked dimensions are bounded using a hyperbolic tangent activation and rescaled to maximum amplitudes of 0\.5 and 0\.1, respectively\. The total log\-likelihood of data sample𝐗¯1,k\\bar\{\\mathbf\{X\}\}\_\{1,k\}, following from Equation \([14](https://arxiv.org/html/2606.30936#S3.E14)\), is then written as

log⁡p​\(𝐗¯1,k\|𝐘k\)=log⁡p​\(𝐰\)\+∑j=1nalog⁡\|𝐉​\(fj\)\|,\\log p\(\\bar\{\\mathbf\{X\}\}\_\{1,k\}\|\\mathbf\{Y\}\_\{k\}\)=\\log p\(\\mathbf\{w\}\)\+\\sum\_\{j=1\}^\{n\_\{a\}\}\\log\\left\|\\mathbf\{J\}\(f\_\{j\}\)\\right\|,\(16\)where the base distribution is chosen as a standard multivariate normalp​\(𝐰\)∼𝒩​\(𝟎,𝐈\)p\(\\mathbf\{w\}\)\\sim\\mathcal\{N\}\(\\mathbf\{0\},\\mathbf\{I\}\), yielding:

log⁡p​\(𝐰\)=−12​𝐰T​𝐰−62​log⁡\(2​π\)\.\\log p\(\\mathbf\{w\}\)=\-\\frac\{1\}\{2\}\\mathbf\{w\}^\{T\}\\mathbf\{w\}\-\\frac\{6\}\{2\}\\log\(2\\pi\)\.

### 3\.2Generative Orbit Determination

The Orbit Determination architecture combines a Conditional Normalising Flow model with a Transformer\-based encoder\[[15](https://arxiv.org/html/2606.30936#bib.bib25)\]\. This class of models is particularly effective for a wide range of tasks, including time\-series analysis, owing to the attention mechanism that enables the extraction of long\-range dependencies within the data\. In the present framework, the Transformer processes the observation sequence

𝐘k\\mathbf\{Y\}\_\{k\}and produces a fixed\-size representation

𝐘~k∈ℝ256\\tilde\{\\mathbf\{Y\}\}\_\{k\}\\in\\mathbb\{R\}^\{256\}\. The purpose of this encoding is to obtain an architecture capable of handling observations originating from different observers and locations, while remaining adaptable to other data modalities \(for example, pre\-processed images\)\. Referring the reader to specialised literature for a detailed description, the Transformer used in this work consists of

nh=16n\_\{h\}=16attention heads and

nl=8n\_\{l\}=8stacked layers\.

Figure[4](https://arxiv.org/html/2606.30936#S3.F4)illustrates the overall Generative Orbit Determination framework\. During training, the model learns the conditional mapping associated with the latent base distribution

p​\(𝐰\)p\(\\mathbf\{w\}\)via the invertible transformations defined in Equation \([12](https://arxiv.org/html/2606.30936#S3.E12)\)\. At inference time, samples are drawn from the latent distribution and mapped through the inverse flow conditioned on new observations, thereby producing estimates of the target state at the beginning of the observation arc\. In this work, sampling is performed in a local neighbourhood of latent codes obtained from the training data, in a manner analogous to a variational setting\. Specifically, perturbed latent samples are generated as

𝐰inf=𝐰train\+0\.1​𝜺,𝜺∼𝒩​\(𝟎6,𝐈6\),𝐰train=f​\(𝐗¯1,k\|𝐘k\),k=1,…,ntrain\.\\mathbf\{w\}\_\{\\mathrm\{inf\}\}=\\mathbf\{w\}\_\{\\mathrm\{train\}\}\+0\.1\\,\\boldsymbol\{\\varepsilon\},\\ \\boldsymbol\{\\varepsilon\}\\sim\\mathcal\{N\}\(\\mathbf\{0\}\_\{6\},\\mathbf\{I\}\_\{6\}\),\\ \\mathbf\{w\}\_\{\\mathrm\{train\}\}=f\(\\bar\{\\mathbf\{X\}\}\_\{1,k\}\|\\mathbf\{Y\}\_\{k\}\),\\ k=1,\\ldots,n\_\{\\mathrm\{train\}\}\.\(17\)Here,𝐰inf∈ℝ6\\mathbf\{w\}\_\{\\mathrm\{inf\}\}\\in\\mathbb\{R\}^\{6\}denotes the latent variable used at inference time, obtained by adding a small Gaussian perturbation to a latent code𝐰train\\mathbf\{w\}\_\{\\mathrm\{train\}\}corresponding to a training sample\. The perturbed latent vectors are then mapped through the inverse transformationg\(⋅\|𝐘new\)g\(\\cdot\|\\mathbf\{Y\}\_\{\\mathrm\{new\}\}\)conditioned on the new observations, yielding a set of plausible initial states consistent with the measurement arc\. A key advantage of this approach is that multiple samples can be generated for the same observation set\. The corresponding outputs constitute a collection of plausible target states consistent with the measurements, enabling uncertainty quantification and improving the robustness of the estimation process\.

![Refer to caption](https://arxiv.org/html/2606.30936v1/Figures/od_model.png)Figure 4:Generative Orbit Determination framework\. Top panel: training process, where dashed lines indicate the backpropagation of the loss function\. Bottom panel: inference stage
### 3\.3Physics\-Informed Loss

The generative model is trained by minimising a loss function defined as the average per sample of the conditional log\-likelihood in Equation \([16](https://arxiv.org/html/2606.30936#S3.E16)\):

ℒflow=−1ntrain​∑k=1ntrainlog⁡p​\(𝐗¯1,k\|𝐘k\)\.\\mathcal\{L\}\_\{\\mathrm\{flow\}\}=\-\\frac\{1\}\{n\_\{\\mathrm\{train\}\}\}\\sum\_\{k=1\}^\{n\_\{\\mathrm\{train\}\}\}\\log p\(\\bar\{\\mathbf\{X\}\}\_\{1,k\}\|\\mathbf\{Y\}\_\{k\}\)\.\(18\)Optimisation of this objective encourages the model to learn an invertible transformation whose induced latent distribution, when mapped back to the state space, reproduces the empirical distribution of initial states conditioned on the observations\. However, likelihood maximisation alone does not guarantee physical consistency\. In particular, the flow may generate initial states that are statistically plausible under the learned distribution but dynamically incompatible with the underlying equations of motion, meaning that time propagation of such states does not reproduce the corresponding observation set\. To enforce dynamical consistency, an additional physics\-informed regularisation term is introduced\. The physics\-informed lossℒphys\\mathcal\{L\}\_\{\\mathrm\{phys\}\}is defined as

ℒphys=1ntrain​∑k=1ntrain‖𝝈¯k‖2,\\mathcal\{L\}\_\{\\mathrm\{phys\}\}=\\frac\{1\}\{n\_\{\\mathrm\{train\}\}\}\\sum\_\{k=1\}^\{n\_\{\\mathrm\{train\}\}\}\\left\\\|\\bar\{\\boldsymbol\{\\sigma\}\}\_\{k\}\\right\\\|^\{2\},\(19\)where𝝈¯k\\bar\{\\boldsymbol\{\\sigma\}\}\_\{k\}denotes the residual vector associated with samplekk, obtained by propagating the candidate initial state forward in time according to the equations of motion and comparing the resulting predicted observations with the measured ones\. The propagation is performed after transforming the initial state to the synodic frame and subsequently back to the geocentric frame in order to compute the topocentric estimated angular observables\. During training, time integration is carried out using a classical fourth\-order Runge–Kutta \(RK4\) scheme, providing a compromise between numerical accuracy and computational efficiency\. The overall training objective combines the probabilistic flow loss and the physics\-informed regularisation through a weighting coefficientβ\\beta:

ℒ=ℒflow\+β​ℒphys\.\\mathcal\{L\}=\\mathcal\{L\}\_\{\\mathrm\{flow\}\}\+\\beta\\,\\mathcal\{L\}\_\{\\mathrm\{phys\}\}\.\(20\)This composite loss promotes solutions that are both statistically consistent with the training data and physically admissible under the dynamical model, thereby encouraging the learned latent manifold to represent meaningful orbital states rather than purely data\-driven interpolations\.

## 4Results

The generative model is trained on the training subset using mini\-batches of up to 128 samples\. Optimisation is performed for 200 epochs using the Adam optimiser with a learning rate of1×10−41\\text\{\\times\}\{10\}^\{\-4\}\. The weighting factorβ\\betalinearly ranges between 0 and 0\.5 in 100 epochs and it is then left constant\. The evolution of the training process is shown in Figure[5](https://arxiv.org/html/2606.30936#S4.F5), where dashed and solid curves correspond to the training and validation losses, respectively\. The blue and red curves represent the two principal components of the flow loss, namely the log\-likelihood of the latent variables under the base distribution and the logarithm of the absolute determinant of the Jacobian of the transformation\.

![Refer to caption](https://arxiv.org/html/2606.30936v1/Figures/Figure_59_1_.png)Figure 5:Training dynamics of the physics\-informed generative modelIt is instructive to compare the learned latent log\-likelihood with the value corresponding to an ideal standard Gaussian distribution in six dimensions:

log⁡𝒩​\(𝐰\)=−62​\(1\+log⁡2​π\)≈−8\.51\.\\log\\mathcal\{N\}\(\\mathbf\{w\}\)=\-\\frac\{6\}\{2\}\\left\(1\+\\log 2\\pi\\right\)\\approx\-8\.51\.The log\-likelihood achieved by the flow model is slightly higher \(less negative\) but remains close to this reference value, indicating that the learned latent distribution approximates a Gaussian while retaining measurable deviations from it\. This behaviour is expected, as the true posterior over physically admissible initial states conditioned on the observations is generally non\-Gaussian\. These deviations justify the use of a local, variational\-style sampling strategy at inference time\. By sampling in the neighbourhood of latent codes obtained from the training data, one effectively exploits informative regions of the latent manifold that correspond to dynamically consistent solutions\. When mapped back through the inverse flow, such samples yield plausible initial states that are both statistically consistent with the learned distribution and physically compatible with the underlying dynamics\.

Table 3:Model training: loss function termsFor ease of reference, Table[3](https://arxiv.org/html/2606.30936#S4.T3)reports a summary of the relevant loss terms for the training and validation sets\.

### 4\.1Estimation Performance Assessment

The primary motivation for adopting a generative approach, as discussed previously, lies in its capability to provide physics\- and context\-informed initial guesses for standard estimation techniques, such as the Nonlinear Least Squares Estimation \(NLSE\) method introduced earlier\. In this section, the performance of the proposed flow\-based model as an initialisation strategy is assessed on the samples of the test subset\. For each observation arc, an initial state estimate is generated following the procedure described above and subsequently used to initialise the NLSE algorithm\. To evaluate its effectiveness, the flow\-based initialisation is compared against a baseline strategy based on the Moon’s state at the corresponding epoch, computed analytically under the same modelling assumptions adopted throughout this work\. This choice is motivated by the proximity of the target orbits to the Moon, which makes such an approximation a plausible practical alternative\. In addition, the exact target state is also used as an initial condition for the estimation process, providing a reference benchmark and ensuring a fair comparison across all cases\. The results are summarised in Figure[6](https://arxiv.org/html/2606.30936#S4.F6)\. The left panel reports the estimation errors in position and velocity \(in normalised units\), while the right panel shows the residuals after convergence of the estimation algorithm\.

![Refer to caption](https://arxiv.org/html/2606.30936v1/Figures/Figure_78.png)Figure 6:Performance assessment of the estimation process initialised with the exact state, the Moon\-based approximation, and the flow model\. Left: position and velocity errors; right: post\-fit residuals\.The results in the left panel indicate that the initialisation provided by the flow model consistently outperforms the Moon\-based approximation and approaches the performance obtained when initialising from the exact state\. This behaviour is observed not only in terms of average performance \(mean values, denoted by dots\) but also in terms of robustness, as evidenced by the reduced spread and the alignment of median values \(squares\) across the test set\. The right panel corroborates these findings\. The residuals obtained using the flow\-based initialisation are of the same order of magnitude as those achieved with the exact initial state, indicating that the estimator converges to physically consistent solutions\. In contrast, the Moon\-based initialisation yields significantly larger residuals and a wider dispersion, suggesting frequent convergence to local minima rather than to the correct solution\.

### 4\.2Physics\-Informed Loss

Inspection of the training dynamics reported in Figure[5](https://arxiv.org/html/2606.30936#S4.F5)might suggest that the physics\-informed loss termℒphys\\mathcal\{L\}\_\{\\mathrm\{phys\}\}plays a marginal role during optimisation, as it remains both small in magnitude and nearly constant throughout the training process\. This behaviour is, to some extent, expected\. Indeed, the normalising flow, by construction, learns an invertible mapping between the data space and the latent space; as a result, it is capable of reconstructing the training samples with high fidelity, implicitly capturing the underlying physical structure of the data throughimitation, even in the absence of an explicit physics constraint\. However, the generative Orbit Determination problem does not merely require reconstruction of known samples, but rather the ability to generate new, physically consistent states when sampling from the latent distribution conditioned on the observations\. This implies that physical admissibility must be enforced not only at the training points, but also in the neighbourhood of the latent manifold where inference is performed\. To further investigate this aspect, the same model was trained under identical conditions but without the physics\-informed loss term\. The corresponding training dynamics are shown in Figure[7](https://arxiv.org/html/2606.30936#S4.F7)\. While the final values of the log\-likelihood components are comparable to those obtained in the physics\-informed case, the absence ofℒphys\\mathcal\{L\}\_\{\\mathrm\{phys\}\}results in a less constrained transformation\. In particular, the physics\-informed loss acts as a regularising mechanism, effectively damping the dynamics of the learned mapping and encouraging the model to preserve the physical consistency of the generated states across the latent space\.

![Refer to caption](https://arxiv.org/html/2606.30936v1/Figures/Figure_81.png)Figure 7:Training dynamics of the physics\-agnostic generative modelTable 4:Model training: loss function termsThe impact of this difference becomes more evident when evaluating the estimation performance obtained using the flow\-based initialisation\. As reported in Table[4](https://arxiv.org/html/2606.30936#S4.T4), although the residuals remain relatively small on average, the estimation exhibits significant degradation in terms of robustness\. In particular, the large mean and maximum errors in both position and velocity indicate the presence of failure cases, in which the estimator converges to physically inconsistent or suboptimal solutions\.

The physics\-informed loss acts as a regularisation term on the latent space geometry\. By enforcing dynamical consistency during training, it ensures that samples drawn from the posterior distribution are not only statistically plausible but also physically admissible\. Consequently, the inclusion of

ℒphys\\mathcal\{L\}\_\{\\mathrm\{phys\}\}is essential to guarantee that the generative model produces reliable initial conditions for downstream estimation tasks, especially in regimes where extrapolation beyond the training data is required\.

## 5Conclusion

This work has demonstrated the application of a generative artificial intelligence framework to the solution of an angles\-only orbit determination problem in the cislunar environment, thereby extending the scope of Generative Astrodynamics\.

Conditional normalising flows have been shown to effectively approximate the probability distribution of the initial state of an observation arc conditioned on the corresponding measurements\. The introduction of a physics\-informed loss, based on the numerical integration of the estimated state and the subsequent application of the measurement model, proved instrumental in regularising the learned transformation\. In particular, it enforces consistency with the underlying dynamical model in the neighbourhood of the latent samples, thus improving both physical fidelity and generalisation\.

At inference time, the proposed approach yields competitive performance when employed as a physics\- and context\-informed initialisation strategy for classical estimation methods such as Nonlinear Least Squares Estimation \(NLSE\)\. In particular, it outperforms initialisation based on the Moon state in terms of both residuals and state estimation errors, while achieving accuracy comparable to that obtained with exact initialisation\. This demonstrates the capability of the method to reliably initialise the estimation process, mitigating the risk of convergence to local minima\.

Moreover, in contrast to deterministic machine learning approaches, the proposed framework estimates a full conditional probability distribution\. This enables the generation of multiple plausible initial states through sampling, supporting more robust estimation pipelines and providing a natural mechanism for uncertainty quantification\. The flexibility of the variational sampling process further allows the representation of complex and potentially multi\-modal distributions, opening the possibility of jointly modelling heterogeneous families of orbits within a unified latent space\.

Owing to its grounding in both physical principles and orbital context, the proposed initialisation strategy has the potential to be extended to a broader class of orbit determination and space navigation problems, including scenarios with sparse or noisy observations, multi\-sensor data fusion, and autonomous onboard navigation\. Future work will focus on assessing the performance of the method in more realistic settings, including real observational data, as well as investigating scalability, robustness to model mismatch, and integration within operational orbit determination frameworks\.

## References

- \[1\]M\. Anand\(2010\-12\)Lunar Water: A Brief Review\.Earth, Moon, and Planets107\(1\),pp\. 65–73\.External Links:ISSN 1573\-0794,[Link](https://doi.org/10.1007/s11038-010-9377-9),[Document](https://dx.doi.org/10.1007/s11038-010-9377-9)Cited by:[§1](https://arxiv.org/html/2606.30936#S1.p2.1)\.
- \[2\]F\. Caldas and C\. Soares\(2024\-07\)Machine Learning in Orbit Estimation: a Survey\.Acta Astronautica220,pp\. 97–107\.Note:arXiv:2207\.08993 \[astro\-ph\]Comment: Accepted for Publication to Acta AstronauticaExternal Links:ISSN 00945765,[Link](http://arxiv.org/abs/2207.08993),[Document](https://dx.doi.org/10.1016/j.actaastro.2024.03.072)Cited by:[§1](https://arxiv.org/html/2606.30936#S1.p4.1)\.
- \[3\]S\. Creech, J\. Guidi, and D\. Elburn\(2022\)Artemis: an overview of nasa’s activities to return humans to the moon\.In2022 IEEE Aerospace Conference \(AERO\),Vol\.,pp\. 1–7\.External Links:[Document](https://dx.doi.org/10.1109/AERO53065.2022.9843277)Cited by:[§1](https://arxiv.org/html/2606.30936#S1.p2.1)\.
- \[4\]L\. Dinh, D\. Krueger, and Y\. Bengio\(2015\)NICE: non\-linear independent components estimation\.InProceedings of the International Conference on Learning Representations \(ICLR\) Workshop,San Diego, CA, USA\.Cited by:[§3\.1](https://arxiv.org/html/2606.30936#S3.SS1.p1.6)\.
- \[5\]L\. Dinh, J\. Sohl\-Dickstein, and S\. Bengio\(2017\-04\)Density estimation using real NVP\.InProceedings of the International Conference on Learning Representations \(ICLR\),Toulon, France\.Cited by:[§3\.1](https://arxiv.org/html/2606.30936#S3.SS1.p2.2)\.
- \[6\]J\. Hartikainen, M\. Seppanen, and S\. Sarkka\(2012\-06\)State\-Space Inference for Non\-Linear Latent Force Models with Application to Satellite Orbit Prediction\.arXiv\.Note:arXiv:1206\.4670 \[cs\]Comment: ICML2012External Links:[Link](http://arxiv.org/abs/1206.4670),[Document](https://dx.doi.org/10.48550/arXiv.1206.4670)Cited by:[§1](https://arxiv.org/html/2606.30936#S1.p4.1)\.
- \[7\]I\. Kobyzev, S\. J\. D\. Prince, and M\. A\. Brubaker\(2021\-11\)Normalizing flows: an introduction and review of current methods\.IEEE Transactions on Pattern Analysis and Machine Intelligence43\(11\),pp\. 3964–3979\.Cited by:[§3\.1](https://arxiv.org/html/2606.30936#S3.SS1.p1.6)\.
- \[8\]W\. Litteri, A\. Francisco Gil, M\. Vasile, V\. Rodriguez\-Fernandez, and D\. Camacho\(2025\-01\)Generative Astrodynamics: Trajectory Analysis and Design in the Restricted Three\-Body Problem\.Cited by:[§1](https://arxiv.org/html/2606.30936#S1.p5.1)\.
- \[9\]W\. Litteri, A\. F\. Gil, M\. Vasile, V\. Rodriguez\-Fernandez, and D\. Camacho\(2026\-01\)Generation of Periodic Orbits in the Restricted Three\-Body Problem with a Variational Autoencoder\.Research Square\.Note:ISSN: 2693\-5015External Links:[Link](https://www.researchsquare.com/article/rs-8651179/v1),[Document](https://dx.doi.org/10.21203/rs.3.rs-8651179/v1)Cited by:[§1](https://arxiv.org/html/2606.30936#S1.p5.1)\.
- \[10\]L\. Manduchi, K\. Pandey, C\. Meister, R\. Bamler, R\. Cotterell, S\. Däubener, S\. Fellenz, A\. Fischer, T\. Gärtner, M\. Kirchler, M\. Kloft, Y\. Li, C\. Lippert, G\. d\. Melo, E\. Nalisnick, B\. Ommer, R\. Ranganath, M\. Rudolph, K\. Ullrich, G\. V\. d\. Broeck, J\. E\. Vogt, Y\. Wang, F\. Wenzel, F\. Wood, S\. Mandt, and V\. Fortuin\(2025\-03\)On the Challenges and Opportunities in Generative AI\.arXiv\.Note:arXiv:2403\.00025 \[cs\] version: 3External Links:[Link](http://arxiv.org/abs/2403.00025),[Document](https://dx.doi.org/10.48550/arXiv.2403.00025)Cited by:[§1](https://arxiv.org/html/2606.30936#S1.p5.1)\.
- \[11\]G\. Papamakarios, E\. Nalisnick, D\. J\. Rezende, S\. Mohamed, and B\. Lakshminarayanan\(2021\)Normalizing flows for probabilistic modeling and inference\.Journal of Machine Learning Research22\(57\),pp\. 1–64\.Cited by:[§3\.1](https://arxiv.org/html/2606.30936#S3.SS1.p1.6)\.
- \[12\]S\. Ross, W\. Koon, M\. Lo, and J\. Marsden\(2022\-10\)Dynamical Systems, the Three\-Body Problem, and Space Mission Design\.External Links:ISBN 978\-0\-615\-24095\-4Cited by:[§2\.1](https://arxiv.org/html/2606.30936#S2.SS1.p1.6)\.
- \[13\]A\. Scorsoglio, A\. D’Ambrosio, L\. Ghilardi, R\. Furfaro, and V\. Reddy\(2023\)Physics\-Informed Orbit Determination for Cislunar Space Applications\.\(en\)\.Cited by:[§1](https://arxiv.org/html/2606.30936#S1.p4.1)\.
- \[14\]D\. Vallado\(1997\-03\)Fundamentals of Astrodynamics and Applications\.External Links:[Link](https://www.semanticscholar.org/paper/Fundamentals-of-Astrodynamics-and-Applications-Vallado/4e2d5e7bdb55431714087aa3160afff8d32b8d5d)Cited by:[§1](https://arxiv.org/html/2606.30936#S1.p1.1)\.
- \[15\]A\. Vaswani, N\. Shazeer, N\. Parmar, J\. Uszkoreit, L\. Jones, A\. N\. Gomez, Ł\. Kaiser, and I\. Polosukhin\(2017\)Attention is all you need\.InAdvances in Neural Information Processing Systems,Vol\.30\.Cited by:[§3\.2](https://arxiv.org/html/2606.30936#S3.SS2.p1.5)\.
- \[16\]S\. Wang, S\. Sankaran, H\. Wang, and P\. Perdikaris\(2023\-08\)An Expert’s Guide to Training Physics\-informed Neural Networks\.arXiv\.Note:arXiv:2308\.08468 \[cs\]Comment: 36 pages, 25 figures, 13 tablesExternal Links:[Link](http://arxiv.org/abs/2308.08468),[Document](https://dx.doi.org/10.48550/arXiv.2308.08468)Cited by:[§1](https://arxiv.org/html/2606.30936#S1.p4.1)\.
- \[17\]J\. Watson, C\. Song, O\. Weeger, T\. Gruner, A\. T\. Le, K\. Pompetzki, A\. Hendawy, O\. Arenz, W\. Trojak, M\. Cranmer, C\. D’Eramo, F\. Bülow, T\. Goyal, J\. Peters, and M\. W\. Hoffman\(2025\-05\)Machine Learning with Physics Knowledge for Prediction: A Survey\.arXiv\(en\)\.Note:arXiv:2408\.09840 \[cs\]Comment: 61 pages, 8 figures, 2 tables\. Accepted at the Transactions of Machine Learning Research \(TMLR\)External Links:[Link](http://arxiv.org/abs/2408.09840),[Document](https://dx.doi.org/10.48550/arXiv.2408.09840)Cited by:[§1](https://arxiv.org/html/2606.30936#S1.p4.1)\.
- \[18\]C\. Wilson and M\. Vasile\(2026\-01\)Generating stable and metastable critical points in uncertain systems via flow‐based models\.Expert Systems43\(3\)\.External Links:ISSN 0266\-4720,[Document](https://dx.doi.org/10.1111/exsy.70196)Cited by:[§1](https://arxiv.org/html/2606.30936#S1.p5.1)\.