Neural network surrogates for elliptic Lambert transfers in preliminary mission design

1 Introduction

Preliminary mission design involves the evaluation of many candidate transfers over launch windows, arrival windows, and candidate gravity assist sequences under coupled timing and propellant constraints (Vasile and De Pascale, 2006). At this stage, the objective is to explore the transfer landscape efficiently and identify feasible regions and low-cost opportunities for subsequent refinement. The models used in this phase must retain sufficient fidelity to capture the main structure of the design space and remain inexpensive enough for repeated use in parameter sweeps and global search procedures. Because Lambert arcs are central to many impulsive transfer analyses (Sangrà and Fantino, 2021), their repeated evaluation can become a limiting factor in early-stage trajectory exploration.

Classical Lambert solvers provide accurate numerical solutions and are a reference approach in impulsive trajectory design (Gooding, 1990; Izzo, 2015). However, in early-stage mission analysis, where very large numbers of candidate transfers must be evaluated across multidimensional search spaces, repeated Lambert evaluations can still represent a considerable computational effort (Wagner et al., 2015).

Recent machine learning approaches in astrodynamics have explored data-driven surrogates and learned initializations for computationally intensive trajectory problems (Izzo et al., 2019). For Lambert-type transfer problems in particular, previous work has mainly employed neural networks either as regression models for transfer parameters or as components of hybrid pipelines in which a learned prediction is followed by numerical refinement (Yang B. et al., 2022).

Beyond the nonlinearity of Lambert’s map, a key difficulty for surrogate learning is that different transfer geometries induce different admissible ranges for the solution variables (Izzo, 2015). The regression problem is not posed over a fixed common domain. The scale and valid interval of the solution variables depend on the specific transfer configuration. For some geometries, the admissible range may be narrow, whereas for others it may be broader and shifted. This makes direct learning more difficult, because the model must account not only for the map itself, but also for the way its valid domain changes from one transfer instance to another. This suggests that improved surrogate performance may depend less on increasing architectural complexity than on constructing a geometry-aware representation that transforms diverse transfer instances into a more regular learning domain.

Motivated by this observation, we study the bounded zero-revolution elliptic Lambert problem under Lagrange’s formulation and examine whether a geometry-aware normalization can regularize the inverse map for surrogate learning. The proposed approach represents each transfer instance in canonical coordinates with bounded, geometry-dependent ranges and uses neural surrogates to predict the semi-major axis of the transfer ellipse. A standard multilayer perceptron is compared with DeepONet and Kolmogorov-Arnold Networks to evaluate predictive accuracy and to examine the structure revealed by the normalized formulation. The aim is to assess whether such surrogates can preserve the main features of the transfer landscape and support the identification of promising candidate regions in preliminary mission design.

2 Lambert's problem and related work

Lambert’s theorem was first formulated by Johann Lambert in 1761 as a problem of finding the orbits of comets based on observational data (Baloglou et al., 2025 citing Lambert 1761). Interest in the theorem was renewed at the start of the space age as a formulation of the problem of finding orbits that intercept a target at a given position and time, which is of high importance when planning rendezvous between two spacecraft or interplanetary transfers. Some notable early works are from Lancaster and Blanchard (1969) and Bate et al. (1971), who both introduced different formulations of Lambert’s equation unifying the cases of elliptic, hyperbolic, and parabolic orbit solutions. Later, Gooding (1990) introduced an efficient method for numerically solving the unified equations presented by Lancaster and Blanchard (1969). Despite its long history, Lambert’s theorem (also known as Lambert’s problem) and methods for efficiently finding solutions remain areas of active research due to their usefulness in orbital dynamics. The classical methods for solving Lambert’s problem can be classified as either iterative or direct methods.

Iterative methods are the most common in the literature and involve two steps: first, an initial guess of the true solution is obtained, and then an iteration scheme is used to iterate on the solution until a sufficiently small error is obtained. Izzo (2015) extended on Gooding (1990) method by introducing auxiliary variables and simplifying the problem, which enabled a more effective way to obtain a good initial guess for the solution to the Lambert problem. By coupling the initial guess with a household iterative scheme for refining the solution, the solver introduced by Izzo (2015) was faster than the previous algorithm by Gooding (1990). Shu et al. (2022) builds on Izzo’s extension to Gooding’s universal variable formulation, by finding higher-order roots using power series-expanded iterative methods. For Earth-Mars transfers, the higher-order roots reduce errors by four orders of magnitude over linearized methods. While speed comparisons are not provided, the approach enables uncertainty analysis for orbit determination and debris cloud evaluation.

Based on Bate et al. (1971) universal variable formulation Arora and Russell (2014) (cited by Russell (2022)) introduced the Vercosine formulation. Russell (2022) then improved the formulation with a new variable transformation that keeps the solution space nearly invariant to the number of revolutions. In the same work, a method for solving the equation is presented, in which initial guesses are generated by interpolating the solution space using multiple basis functions (requiring precomputed coefficients), and a root solver refines the solution. The method achieves accuracy across the full solution space, with reported speedups exceeding 2 × over Gooding’s method, while adding the requirement of storing the precomputed coefficients.

De La Torre et al. (2018) presents a reformulation of Lambert’s problem using Levi-Civita regularization. The formulation is extended to handle multi-revolution orbits, and a solution-finding method is presented that uses lookup tables for the initial guess, which are then refined through Newton-Raphson iterations. Recent work by Pontani et al. (2026) formulates the problem in more traditional Keplerian terms, focusing on finding the true anomaly for multi-revolution elliptic orbits.

Direct methods are less common but focus on finding satisfactory solutions without the iterative step. Bombardelli et al. (2018) presented an analytical method for finding approximate solutions to Lambert’s problem. Negrete and Abdelkhalik (2025) introduced a method based on Bate’s universal variable formulation and complex analysis. It applies contour integrals near the true solution, with the accuracy of the results depending on how well the integrals can be computed. The reported execution time is faster than previous approaches, including Russell (2022) and Gooding (1990).

In recent years, deep learning-based approaches have started to gain traction in astrodynamics research, where neural networks have been used to find low-thrust trajectories in the asteroid belt (Acciarini et al., 2024), designing flyby asteroid sequences (Ozaki et al., 2022), or computing flyby maps in the Jupiter moon system (Yang H. et al., 2022), to name a few examples.

Regarding Lambert’s problem, Choi et al. (2024) evaluates the use of a basic Multilayer Perceptron (MLP) network to approximate solutions to Lambert’s Problem, approximating the velocity of the transfer orbit at the initial position. In more recent work, Choi et al. (2025) presents a hybrid method: first, an MLP approximates the universal variable as in the procedure given by Gooding (1990); if higher precision is needed, the solution is refined using Halley’s method. Gueho et al. (2020) compared parametric and machine learning-based approaches to study the uncertain Lambert problem. Solving Lambert’s problem yields a purely Keplerian orbit; however, in practice, additional perturbations, such as the J2 zonal harmonics arising from the central body’s non-uniform gravity field, must also be considered. Yang B. et al. (2022) proposed a neural network to enhance the Lambert solution for the J2-perturbed problem in the Jupiter system by providing an initial guess for the J2-perturbed case, accelerating convergence to the final solution.

3 Governing equations

Lambert’s problem is a boundary value problem where you want to find an orbit between two positions, with a specific time-of-flight between them. In this work, we will use a formulation of Lambert’s problem based on Lagrange’s solution (Baloglou et al., 2025). For elliptical orbits, this can be expressed with the equation,where t t − t i is the time difference between the initial and target positions, i.e., the time-of-flight ( τ ) as we will refer to it hereafter. In this equation, we have one constant, μ - the standard gravitational parameter, and two unknowns, a - the semi-major axis of the ellipse, and Q - the number of revolutions around the orbit before reaching the target position. There are also two auxiliary variables, α and β, which are defined in terms of the position vectors (r i ⃗, r t ⃗) and the semi-major axis. where c is the chord and s is the semiperimeter of the triangle formed by the ellipse’s focus point corresponding to the parent body and the positions r i and r t.

t t − t i = a 3 μ 2 π Q + α − sin α − β − sin β (1)

α = α i , longway solution 2 π − α i , shortway solution (2)

α i = 2 ⁡ sin − 1 s 2 a

β = 2 ⁡ sin − 1 s − c 2 a

c = ‖ r t ⃗ − r i ⃗ ‖

s = c + ‖ r i ⃗ ‖ + ‖ r t ⃗ ‖ 2 (3)

As can be seen from these equations, and in Figure 1b, there are multiple τ for a given a, first depending on how many orbits around the parent body we want, Q in Equation 1, before reaching the target location ¹ . However, also depending on what is called the shortway or the longway solution, Equation 2. The shortway or longways solutions stem from the fact that there exist two different ellipses, with the same semi-major axis a and focus point F, which intersect the two positions r i and r t, see Figure 1a. From Figure 1b, we also see that there is a minimum value of the semi-major axis, a min. This is the smallest ellipse that connects the two positions and can be found with

a min = s 2 (4)

Solving Lambert’s problem involves the inverse of Equation 1, as we want to find a given r i, r t, and τ. However, Equation 1 is not analytically invertible and numerical methods, such as the one presented by Izzo (2015), are required to find the semi-major axis. In practice, this means that we can only find an approximate semi-major axis ( a ̃ ) from our numerical inversion F n − 1:

a ̃ = F n − 1 τ , r 1 ⃗ , r 2 ⃗ (5)

3.1 Multi-planetary flyby

Multi-planetary flybys are used in interplanetary mission design to satisfy low delta-v budgets and reduce propellant requirements by leveraging gravitational assists, allowing trajectories to reach more distant targets. Such trajectories are generally modeled as a chain of heliocentric transfer arcs, each obtained by repeatedly solving Lambert’s problem between consecutive planetary encounters. When assessing multi-planetary flyby maneuvers, the incoming and outgoing velocities in the planetary reference frame for each flyby can be computed as follows where v + and v − are the spacecraft velocity vectors in the solar reference frame, before and after the flyby, respectively, and v p is the planet’s velocity in the solar reference frame. The incoming velocity, v −, can be obtained from the Lambert solution intercepting the planet; however, v + is unknown and depends on the geometry of the flyby. By assuming a non-powered gravity assist, we know that the size of the incoming velocity in the planet’s reference frame has to be equal to the size of the outgoing velocity (Curtis, 2021)

v ∞ + = v + − v p (6)

v ∞ − = v − − v p (7)

‖ v ∞ + ‖ = ‖ v ∞ − ‖ = v ∞ (8)

If we also assume that after the planetary flyby, the velocity vector of the spacecraft should point along the velocity vector v i − corresponding to the direction of the velocity vector for the transfer to the next planet in the flyby sequence. We can calculate the direction of this vector as the unit vector

v ̂ ∞ − = v i − − v p ‖ v i − − v p ‖ (9)

and can therefore calculate the desired velocity vector according to

v ∞ − = v ∞ v ̂ ∞ − (10)

With Equations 6-10, we can calculate the desired flyby geometry and the closest approach to the planet, thereby determining the outgoing velocity vector. where δ is the turning angle, e is the eccentricity of the hyperbolic orbit of the flyby, and r p is the periapsis, i.e., closest approach (Curtis, 2021). To ensure that your spacecraft does not crash into the planet, r p should be larger than the planetary radius ( R p ) and outside most of the planetary atmosphere, if one exists.

δ = arccos v ∞ + ⋅ v ∞ − v ∞ 2 (11)

e = 1 sin δ / 2 (12)

r p = μ p v ∞ 2 e − 1 (13)

4 Methodology

This section describes the proposed surrogate-learning framework for the bounded zero-revolution elliptic Lambert problem. We first recast the problem in canonical coordinates and construct a dataset over the bounded elliptic solution domain. We then introduce the geometry-aware normalization that maps transfer instances with geometry-dependent admissible ranges into a common representation. Finally, we present the neural surrogate models and the training and evaluation setup used to assess their performance.

4.1 Canonical reduction of the elliptic Lambert problem

To make the bounded elliptic Lambert problem more amenable to surrogate learning, we first express it in canonical units and reduce the three-dimensional geometry to a two-dimensional representation in polar coordinates. Examining Equation 1, we see that it contains the standard gravitational parameter μ, which is dependent on the parent body being orbited. The solution to Equation 1 therefore depends on the parent body. To remove this dependence on the parent body, we can convert the physical units into canonical units. In astrodynamics, this is done by defining a reference distance ( R r ) and setting μ = 1. For example, when working with orbits around Earth, the reference distance is often set to the Earth’s radius (R r = 6378 km). When working with orbits around the Sun, the astronomical unit is often used (R r = 1 ( AU ) = 149,597,870.7 (km)). In this work, we set the reference distance to be equal to the Euclidean norm of r i ⃗. With this, distances, time, and velocities can be converted from physical units to canonical units using the equations. where μ p is the physical μ for the parent body of the system, D, T, and V are the distance, time, and velocity in our original reference frame.

D c = D R r , R r > 0 (14)

T c = T ⋅ μ p R r 3 , R r > 0 , μ p > 0 (15)

V c = V ⋅ R r μ p , R r > 0 , μ p > 0 (16)

We can further simplify the 3D Lambert’s problem into a 2D representation by defining a plane using the r i, r t, and the location of the parent body, which is the focal point of the orbital ellipse (Choi et al., 2024). We convert r t to polar-coordinates in a plane by defining r i to be ( r , θ ) = ( ‖ r i ⃗ ‖ , 0 ). Due to our usage of canonical units relative to r i, we have that ‖ r i ⃗ ‖ = 1. The network approximation for the inverse of the Lambert equation, Equation 5, can therefore be expressed aswhere r and θ describe the target position in polar coordinates with the origin at the orbital ellipse focus and the initial position at r = 1 and θ = 0.

a ̃ = N r , θ , τ (17)

4.2 Data generation for the bounded zero-revolution elliptic domain

Because the inverse Lambert map is not available in closed form, we construct the dataset by sampling the bounded elliptic solution space in the forward direction, using the semi-major axis as the generative variable and computing the corresponding time of flight. While finding the solution to Equation 1 given two positions and a time-of-flight ( τ ) requires numerical methods, finding the τ from the two positions and a semi-major axis ( a ≥ a min ) is trivial. We also note that Lambert’s problem can be formulated as a two-dimensional problem relative to the unit vector r i = ( 1,0 ), as described earlier. With this method, we first generate P target vectors r t ⃗ in polar coordinates relative to r i ⃗, we then generate N semi-major axis values uniformly distributed between a min and a max, a min is defined by Equation 4 and a max is chosen to be equal to s from Equation 3. These values are then used to calculate the τ shortway and the longway solution using a 2D polar version of Equation 1; thereby, there are a total of 2 N − 1 samples for each target vector, we subtract one to account for the fact that the longway and shortway solution for a min are the same.

The generation of r t values is performed by uniform sampling on the r-range [0.1,10] and θ-range [0.1,2 π − 0.1] and all calculations are done using μ = 1. When θ approaches 0 (or 2 π), a in Equation 1 approaches infinity. To prevent this, we restrict θ such that it remains at least 0.1 away from these boundary values. By generating data for r and θ directly, we create a dataset in the canonical units described previously, Equations 14-16, and no conversion is required during network training.

The generated data is then split into a training set and a validation set with a [0.7,0.3] split. Due to the shape of the input and output of the DeepONet architecture, using all of the τ-values for a specific target to give the complete solution for that target, the split is performed based on the target vectors; i.e., all samples with the same r t are assigned to the same set. The parameters used to generate the data are summarized in Table 1.

TABLE 1. Data generation parameters, s is defined in Equation 3.

We evaluate the distribution of the data used for training by plotting the distribution of data points in the training set for different target coordinates (r and θ), the semi-major axis ( a ), and the time-of-flight ( τ ), as shown in Figure 2. As r and θ were sampled uniformly at random from a given range, these distributions (top-left and top-right in Figure 2) are uniform. There are variations across bins, primarily due to the relatively small sample size ( P = 1000 ). The semi-major axis is defined on a uniform distribution; however, as the range for this distribution is defined by the target coordinates, the distribution of the semi-major axis (bottom-left in Figure 2) is not entirely uniform. Between a semi-major axis of 1 to about 5 CDU, multiple target positions have overlapping ranges for the semi-major axis, and the distribution is mostly uniform. For semi-major axis values greater than 5 CDU, fewer targets have overlapping ranges (see Figure 3), and the distribution decreases. A similar effect is present for semi-major axis values smaller than 1. However, for targets with small target radial coordinates, the smaller range between a min and a max leads to tighter sampling; see the left plot in Figure 3. This tighter sampling is responsible for the peak in the data distribution around a semi-major axis of 1 CDU. As with the semi-major axis, the time-of-flight distribution is affected by overlapping ranges and tighter sampling, resulting in a higher concentration of data points for shorter times of flight.

4.3 Geometry-aware normalization of the Lambert map

A challenge in learning the bounded elliptic Lambert map is that the admissible ranges of time of flight τ and semi-major axis a vary with the transfer geometry. To address this challenge, we introduce a geometry-aware normalization that rescales each transfer instance using geometry-dependent bounds derived from the Lambert formulation itself. As illustrated in Figure 3, these geometry-dependent ranges for τ and a vary considerably between different target positions; for example, the valid τ and a for A are [6, 70] and [2.5, 5.3], while C has τ [24,145] and a [4.1, 8.2]. Such variability presents additional challenges for neural networks, especially for architectures like DeepONet, which are designed to learn functional relationships over fixed input domains. In the classical Lambert formulation (Equation 1), a is only bounded below by the minimum a min (Equation 4), and τ is limited by t min as a approaches infinity. However, in our data generation, we have also introduced an upper limit on the semi-major axis ( a max ) based on the generated data range; see Table 1.

Based on a max, a t min for our generated data can be calculated as the shortway solution to Equation 1 given the inputs r i, r t, and a max, where r i and r t are specific to the sample. For our generated data, we use a version of Equation 1 assuming a fixed r i, and using r and θ as input. Using the a min value, we can also calculate the τ for the boundary point between the shortway and longway solutions; here, we will call this t mid. Using these values, calculated for each sample from the input data, we can normalize the input τ and a according to the equations effectively transforming τ = t min → − 1, τ = t mid → 0 and a norm to the range [0,1], see the normalized data in the right plot of Figure 3.

τ norm = τ − t min t mid − t min − 1 (18)

a norm = a − a min a max − a min (19)

To investigate how this geometry-aware normalization affects the training, MLPs with 6 layers, each with 30 nodes and SiLU activation functions, are trained using four different normalization techniques. The first two are standard z-score and minmax normalizations, computed on the training set. The second two are geometry-aware methods, geo-minmax and geo-midmin. geo-minmax applies minmax scaling to a based on Equation 19, and to τ using short- and long-way values for a max. geo-midmin is the normalization outlined earlier which uses Equations 18, 19 to further exploit the underlying geometry of the problem. Figure 4 shows the loss when training the network for 1000 epochs with the ADAM optimizer followed by 400 epochs with the LBFGS optimizer. The geometry-aware normalizations improve the final loss by two orders of magnitude compared to z-score normalization and one order of magnitude compared to minmax.

Leveraging the problem’s geometry in the normalization preserves essential information while making it more suitable for neural networks. As demonstrated in Figure 3, this approach brings the solutions at different positions r i and r t into closer proximity and reveals patterns in the data. As will be demonstrated in Section 5 these are patterns that the networks can learn. The normalized data vary very little, particularly for τ < 0, see Figure 3. However, there is some variation in the time-of-flight for semi-major axis of 1; we refer to these as the τ max values, as they are the maximum τ for a given target. To understand how the data vary across different targets, we can plot the τ max values (x-axis) against the target coordinates, r (y-axis), and θ (colorbar). In Figure 5, we show the distribution of normalized τ max values for the training and validation sets. We observe that as r approaches 10, the τ max values converge to the central value ( ∼ 8.5 ). There is also a similar convergence as r approaches zero. The most spread in the τ max values are around r ≈ 1, however only high (approaching 2 π) or low (approaching 0) target angles θ. For target angles θ ≈ π, the values are still along the central value.

4.4 Neural network surrogate models

Once the bounded Lambert map is expressed in normalized canonical coordinates, we evaluate whether lightweight neural surrogates can learn the resulting inverse relation expressed by Equation 17. We compare a standard multilayer perceptron with two more structured architectures, DeepONet (Lu et al., 2021) and KAN (Liu et al., 2025), in order to assess predictive performance on the normalized problem and to examine what regularities are exposed by the proposed representation.

4.4.1 Multilayer perceptrons

The MLP serves as the primary predictive baseline in this study. MLP networks have long been the workhorse in deep learning. First introduced by Rosenblatt (1958), MLPs have been successfully used in both classification and regression. The principles of the MLP have also been adopted in more advanced neural network architectures, where they operate in tandem with other components. MLPs can be described as a series of linear transformations ( Z l ) and non-linear activations a l transforming the input x to the output y:where l is the network layer (Goodfellow et al., 2016). The linear transformation Z l consists of a weight matrix W l and a bias vector b l, which are applied to the layer input according to

y = a l ◦ Z l ◦ a l ◦ Z l ◦ … a l ◦ Z l ◦ x

y l = Z l ◦ x l = W l x l + b l

in our case the output y will be the approximation of the semi-major axis, a ̃, and the input x is a vector with the target positions and the time-of-flight.

4.4.2 Deep Neural Operator Networks (DeepONet)

Deep Neural Operator Networks were first proposed by Lu et al. (2021). This architecture splits the network into a “Trunk net” and one (or more) “Branch Net(s)”, as shown in Figure 6.

The main idea behind the DeepONet architecture is to use a neural network to approximate operators mapping from one function space to another (Lu et al., 2021). The operator approximation is performed by feeding the input functions into a neural network, the branch net, which estimates the operator coefficients corresponding to the output function at a specific location. In our case, the input to the branch net will be the target position. The DeepONet architecture also introduces the trunk net, which maps the location variable, in our case τ, to the operator coefficients, yielding a set of characteristic equations for the solution. The outputs from the branch and trunk nets are then combined according to

N r i , r t , τ = G r i , r t τ = ∑ k = 0 P b k r i , r t t k τ + b 0

Where b k and t k are the outputs from the branch and trunk nets, respectively, and b 0 is a bias term added to reduce the generalization error (Lu et al., 2021).

4.4.3 Kolmogorov–Arnold networks (KAN)

KAN is the third surrogate architecture considered in this work. The KAN architecture differs significantly from previous architectures in that it is not based on one or more MLP networks and does not rely on the universal approximation theorem; see Figure 6. Instead, the KAN introduces a network based on the Kolmogorov-Arnold Representation theorem, where the network nodes have learnable activation functions instead of learnable input weights (Liu et al., 2025). The KAN achieves this by introducing an activation function on each edge between nodes, see Figure 6, where each of these functions is represented as a sum of B-spline functions ϕ i, multiplied by a coefficient ( c i ):similar to the grid in the finite element method. The coefficients c i are learnable parameters that enable the network to learn the activation functions. This construction increases functional flexibility and provides an alternative structured parameterization of the normalized inverse Lambert map.

Φ x l = ∑ i = 0 S c i ϕ i x l

4.5 Training and evaluation setup

The models are trained on the normalized bounded elliptic dataset and evaluated both on the validation set and in downstream preliminary mission-design scenarios. For the MLP network, we perform hyperparameter tuning by training different models with varying numbers of layers, layer widths, and activation functions for 200 epochs using the Adam optimizer. The various parameters tested are summarized in Table 2. The best configuration was selected and trained for 1000 epochs with the Adam optimizer, followed by 300 epochs with the LBFGS optimizer.

TABLE 2. Hyperparameters evaluated for MLP-based networks.

To create the DeepONet architectures, the best MLP network configuration was used, then split evenly across the subnetworks. Then the networks were trained for 1000 epochs with the Adam optimizer and 300 epochs with the LBFGS optimizer. For DeepONet, we also evaluated different numbers of outputs from the branch and trunk nets, corresponding to varying numbers of characteristic functions generated by the trunk net.

The methodology for training and evaluation of the KAN differs significantly from the MLP-based networks. Because KANs consist of learnable activation functions, they are inherently geared toward function discovery. In this work, we have used the PyKAN framework ² developed by Liu et al. (2025). The PyKAN framework offers a suite of techniques for refining model architecture, not only by increasing complexity through larger activation-function grids or by expanding the network’s depth and width, but also by enforcing sparsity and applying pruning mechanisms. These capabilities support the construction of more parsimonious models, either through a bottom-up approach, by starting with a simple model and adding complexity where it is necessary, or through a top-down approach, by starting from a large model and eliminating connections whose contributions to the outcome are minimal. Once a model has been found, PyKAN supports fitting the activation functions to simple symbolic functions, such as polynomials, sines, or exponentials, enabling the discovery of symbolic formulations of the model (Liu et al., 2024).

For our problem, the top-down approach frequently produced overly complex models. Because our problem is sensitive to slight variations in the final output, we cannot rely on strong regularization to stabilize the solution. However, the top-down method depends heavily on strong regularization of connections. We have therefore used a bottom-up approach in this work, starting from the simplest possible KAN model (see Figure 7), with a single learnable univariate function for each input, followed by a sum of the outputs of the three functions. To preserve a high degree of simplicity in each edge function, we use a second-order B-spline representation with a grid resolution of six. In our experience, this configuration generally yields edge functions that remain easy to interpret and are fast to train, while still providing sufficient flexibility and expressiveness to model the required behavior.

Figure 7 shows an example of the simplest KAN model trained on the data. We can already see some basic structure corresponding to the structures in the data. The time-of-flight, τ, has the basic shape from the normalized data, Figure 7, and the outline of the θ edge function indicates the importance of small and large values, i.e., close to 0 and 2 π. However, r and θ have a lower impact on the final output, as indicated by the faded edges, and the model primarily captures the data’s mean behavior. To capture more structure, we can first expand the model in width, thereby, in analogy with the Kolmogorov-Arnold representation theorem, allowing the model to capture finer-grained structure in the data. We can also expand the model in depth to enable more advanced combinations of our single-variate functions.

4.5.1 Mission-design case studies

In addition to validation-set performance, the trained surrogates are evaluated in two mission-design case studies representative of preliminary trajectory analysis. To evaluate the network on real use-cases, we investigate an Earth-Mars transfer window during 2018, finding the date for the transfer with the lowest delta-V requirement, and a multi-planetary flyby, calculating the delta-V requirement for interplanetary transfers based on the pre-optimized trajectory reported by Gad and Abdelkhalik (2011). We have used a pre-optimization of the system tour without deep-space maneuvers, see Table 3. The network is used to calculate the trajectory arc between each step in the system tour. The required Δ v after each flyby is then calculated as described in Section 3.1 with the assumption that the outgoing velocity vector aligns with the required velocity vector, thereby only requiring an increase in velocity. To obtain the states of the planetary bodies at given times, we used SpicePy (Annex et al., 2020) with Spice Kernels obtained from NASA’s Navigation and Ancillary Information Facility (Acton, 1996).

TABLE 3. Flyby sequence used as test case, retrieved from Gad and Abdelkhalik (2011).

5 Results

This section first identifies one representative configuration for each surrogate architecture, then compares their predictive accuracy on the normalized inverse Lambert map, and finally evaluates their behavior in mission-design case studies.

5.1 Selection of surrogate models

We first identify one representative configuration for each surrogate architecture in order to compare their performance on the normalized inverse Lambert map.

5.1.1 MLP hyperparameter search

We perform a hyperparameter evaluation to identify the best MLP configurations by training the networks as described in Section 4.5 and evaluating the MSE loss on the validation set. Selecting configurations with the lowest MSE, we found that three-layer MLPs with 50 nodes per layer and SiLU activation performed best. Second best used six layers of 30 nodes with SiLU activation functions. Of the ten best network configurations, only the 10th used an activation function other than SiLU, in this case, the Tanh function. However, the MSE loss on the validation set for the 10th best network is one order of magnitude larger than that of the best network configuration. Adding a Sigmoid activation to the network output is an attractive option, as it would produce outputs within the valid range of our normalization. However, our evaluation found that this resulted in higher training loss. Based on this, we selected the two best-performing MLP networks for complete training with the Adam and LBFGS optimizers.

5.1.2 DeepONet model selection

To create the DeepONet, we base the branch and trunk subnetworks on the selected MLP configurations and vary the number of characteristic functions to assess the operator decomposition required by the normalized Lambert map. The top two configurations from the MLP network were selected and split into two equally sized networks. Additionally, because most of the deviation between solutions stems from differences in the targets, we also train networks in which the branch net is larger than the trunk net. In these cases, we split the MLP configurations unevenly, favoring the branch net.

Evaluating the top-5 DeepONet configurations with the lowest MSE loss on the validation set, we find that all use six hidden layers, each with 15 nodes for both the branch and trunk nets. The difference between these networks lies in the number of characteristic functions. Among these, the DeepONet with 10 characteristic functions achieves the lowest loss, but most networks with varying numbers of characteristic functions exhibit similar final losses. For further evaluation, we selected the configuration with two networks, each with six hidden layers of 15 nodes and five characteristic functions. This configuration has the fourth-lowest loss; however, because it uses only five characteristic functions, it is more interpretable than a larger number of functions.

For evaluation of the DeepONet output and how the different networks build the final prediction, we plot the trunk net output (left), the z-score normalized branch net output (center), and the combination (right) in Figure 8. For the trunk net output, we observe that the functions 0, 1, and 3 are very similar to the normalized data in Figure 3. The main difference is that Trunk Net 0 is inverted; however, the branch net output corresponding to Trunk Net 0 has a negative mean, indicating that the Trunk Net 0 output will be inverted in the final DeepONet output. This inversion can be seen in the right-most plot in Figure 8. The outputs for Trunk Net 2 and 4 are very different; firstly, both have a discontinuity in the normalized τ range [ − 1,0 ]. Secondly, Trunk Net 2 has a decreasing slope for normalized τ values greater than 5.

We evaluate how the trunk nets combine to form the final output by assessing the branch net outputs. In the central plot of Figure 8, we show how the branch net outputs depend on target radius r (y-axis) and target angle θ (colorbar). The values are z-score normalized for comparison, with the mean and standard deviation listed in the legend. We observe here that the branch net outputs follow a plot similar to that of the normalized τ max in Figure 5, with the most significant deviation from the central line occurring at around r = 1. Examining the mean and standard deviation values, we observe that the branch net outputs 0, 1, and 3, corresponding to the same-numbered trunk nets, have large means with smaller standard deviations. The branch net outputs corresponding to the trunk nets 2 and 4 have small means and large standard deviations.

5.1.3 KAN model construction

For KAN, a bottom-up construction strategy was adopted, starting from the simplest admissible architecture (Figure 7) and increasing width or depth only when necessary to capture additional structure in the normalized data. We found that, when using the simplest one-layer KAN model with a single edge function per input, the time-of-flight τ dominates the output, and the output varies little with the target positions r and θ. By increasing the network’s depth, we enable more complex combinations of inputs; by increasing its width, we allow a wider range of input combinations. We observed that, with a two-layer network in which the first layer has three nodes and the second layer has a single node (see the left panel of Figure 9), the KAN model achieves sufficiently good performance to yield usable results.

KAN enables fitting the learned functions to symbolic expressions, facilitating function discovery. However, as demonstrated in the right panel of Figure 9, solutions derived from these symbolic formulas deviate substantially from those produced by the KAN model itself, leading to a marked decrease in accuracy. Consequently, the symbolic model is not suitable for solving Lambert’s problem. A more detailed evaluation can be found in the Supplementary Material.

5.2 Comparison of the selected surrogate models

After selecting one representative configuration for each architecture (Table 4), we compare their performance on the validation set using denormalized error in the semi-major axis. The L1 error is computed as a percentage of the actual value and the error is plotted for each model, Figure 10, color-coded according to r (left) and θ (right). The first observation we make is that the MLP has the lowest overall error, followed by DeepONet and then KAN.

TABLE 4. Network configurations selected for final evaluation; see percentage of l1 error on validation set for the trained models in Figure 10.

Evaluating the error for different r values we can see that the higher errors are at smaller r values, the left plots of Figure 10. The areas with smaller r values also correspond to the area of higher variance in solutions observed in Figure 5. We also observe that the percentage error tends to be smaller for larger τ and r. However, it is essential to note that for larger τ and r, the semi-major axis of the transfer orbit is also larger, resulting in a larger absolute error for large r.

Evaluating the error for the different θ values we do not find a clear dependence between the error and the target angle θ for the DeepONet, right plots in Figure 10. However, for the KAN network, it is clear that it has difficulty with θ values near 2 π or 0, i.e., at the extremes where the solution deviates most. We observe a slight tendency toward the same pattern of difficulties for extreme θ values in the MLP network. However, the MLP has significantly less severe errors than those for the KAN.

5.3 Mission-design case studies

To assess whether the selected surrogates remain useful beyond validation-set metrics, we evaluate them in two mission-design case studies representative of preliminary trajectory analysis. The calculation of an Earth-Mars transfer, creating a porkchop plot over the Delta-V requirements for the 2018 transfer window with the identification of the date with the smallest Delta-V requirement. We also evaluate the use of the networks to calculate the Delta-V requirement for a multi-gravity-assist trajectory.

5.3.1 Earth-Mars transfer window

The first case study considers an Earth–Mars transfer window during 2018, where the surrogate-predicted semi-major axis is used to reconstruct transfer velocities and generate porkchop plots for comparison with a classical Lambert solution.

The overall structure of the various porkchop plots in Figure 11 are the same. The most notable differences among the plots occur in the time-of-flight range between 200 and 300 days, which is particularly evident in the DeepONet and KAN plots on the right. In this region, a discontinuity appears in the ML model solutions, corresponding to cases where the semi-major axis approaches its minimum value ( a min ). This discontinuity arises primarily from the post-processing step, in which transfer orbit velocities are computed from positions and semi-major axes.

The logic for calculating the orbital velocities does not use τ in the actual calculation. It does, however, use τ to know if the solution is a shortway or a longway solution. The use of τ leads the velocity calculations to be performed with the assumption of a short time of flight, i.e., Δ τ < 0, until the time-of-flight passes the break point between the shortway and the longway solutions, τ > τ mid, after which the calculations are instead performed for a longer time-of-flight, Δ τ > 0, see Figure 12.

The error is more pronounced on the side for the shortway solution, τ < τ mid, because close to a min on this side, i.e., τ < τ mid, a small error in the estimated semi-major axis ( a est ) give a larger error in the corresponding time-of-flight ( τ est ), than it would on the longway side.

For the MLP model, the region around a min has also yielded numerous invalid solutions, as evidenced by the gray speckles clustered around the discontinuity in the DeepONet and KAN plots. Here, the MLP model returns values for the semi-major axis below a min, resulting in invalid velocities during calculation.

We estimate the level of error from the network methods by calculating the absolute difference between the network solutions and the classical solution, see Figure 13. While the classical solution method is not exact, it has a generally higher accuracy (Izzo, 2015) than have been displayed by the network-based methods in Figure 10. While the error of the MLP network is generally below 1 0 − 2 km/s, it can reach as high as 1 km/s near the a min solution, in the region where we also find invalid solutions. The DeepONet solution has a similar error level, while the highest error, around a min, reaches up to 3.5 km/s. Analogous to Figure 10, the KAN network has the highest error of the three networks. However, the largest error is actually slightly less than for DeepONet. Figure 13 indicates that all of the network models struggle to give accurate solutions around a min.

To identify the date with the minimum delta-V requirement, we computed the minimum value in the delta-V matrix for each method. Table 5 lists the solutions for each of the classical solutions and each of the models. Here, we observe that the MLP model identified the exact date and a delta-V requirement within ∼ 3m/s. The DeepONet and KAN models did not perform as well and estimated earlier departure and arrival dates, with a ∼ 22m/s and ∼ 40m/s higher delta-V requirement, respectively.

TABLE 5. The Departure date, Arrival date, and Delta-V for minimum Delta-V Earth-Mars transfer with departure during 2018 as identified by a classical (Izzo et al., 2020), MLP, DeepONet and KAN.

5.3.1.1 Computational performance

Although comprehensive performance benchmarking was not conducted, a preliminary comparison was carried out with the Lambert solver described by Russell (2022) for the Earth–Mars transfer scenario. The measurements were performed on an Intel i9-10920x CPU at 3.50 GHz with 24 virtual cores and a NVIDIA GeForce RTX 2060 GPU with 8 GB RAM. The solver from Russell (2022) is implemented in Fortran with built-in batching and multi-threading via OpenMP, enabling efficient processing of large trajectory batches. However, because the solver lacks a Python API, a custom Python wrapper was developed to enable execution and benchmarking in Python. As expected, for a small number of trajectories (e.g., a 100 × 100 grid), the traditional solver outperforms the neural network methods. However, when evaluating large-scale problems such as a 3000 × 3000 grid, the GPU-based MLP network achieves a modest speedup of 1.15 ( ± 0.008 ) over the traditional solver, with DeepONet and KAN networks showing slightly greater speedups. Given that neural networks also exhibit higher error, we measured the execution time of the traditional solver with iterative refinement disabled, using only the initial approximation. In this scenario, the network speedup dropped to 1.10 ( ± 0.008 ). It is important to note that only the computation phase was measured, after all r i, r t, and τ values had been assembled in arrays. Additionally, the time measurements for both traditional and neural network methods included the necessary transformation from physical units to each method’s required units. The Russell (2022) solver would likely perform better if the preprocessing was entirely implemented in Fortran, but the same holds for the neural network methods. More recent work (Negrete and Abdelkhalik, 2025) has reported speedups over the method of Russell (2022).

5.3.2 Multi-planetary flyby test case

The second case study evaluates the selected surrogates on a multi-planetary flyby sequence from Earth to Saturn, using fixed encounter dates and neglecting deep-space maneuvers. When planning interplanetary trajectories for uncrewed spacecraft, it is common to use planets and other bodies for gravity assists. Gravity assists can reduce the delta-v requirements for spacecraft, especially when sending them to the outer solar system (Gad and Abdelkhalik, 2011). To evaluate the models for this scenario, we used a trajectory from Earth to Saturn that passes by Venus, Earth, and Jupiter, without deep-space maneuvers; see Table 3 for the full sequence with flyby dates.

The closest approach (periapsis) and delta-v requirements for each flyby is calculated according to Equations 6-13 in Section 3.1 and listed in Table 6 for the classical, MLP, DeepONet, and KAN methods. The total delta-V cost estimated by the ML models is very close to the classical solution, with errors of ∼ 25m/s for the MLP model, ∼ 8m/s for DeepONet, and ∼ 12m/s for KAN. In this example, DeepONet outperforms the MLP. However, when we evaluate the individual transfer and flyby calculations, we observe that the DeepONet model tends to yield higher errors. However, because it sometimes overestimates and sometimes underestimates, the final error is lower for the DeepONet model.

TABLE 6. Comparison of the solar system tour, calculated using a classical (Izzo et al., 2020), MLP, DeepONet, and KAN method.

Here, we also observe the effect of the scales on the absolute error. For example, the Jupiter flyby periapsis calculated using DeepONet has a relative error of only 0.8%. However, in absolute terms, this corresponds to an error of 41,077 km. These results indicate that the surrogates can preserve coarse sequence-level trends, but that individual trajectory elements may still exhibit substantial error.

6 Discussion

The main result of this study is that the bounded zero-revolution elliptic Lambert map becomes substantially more amenable to surrogate learning when expressed through a geometry-aware normalization. The solutions to Lambert’s problem are inherently complex, with edge cases in which the solution tends to infinity, such as when the angle between the initial and target positions approaches zero. However, by imposing limits on the problem, we can create structure to the data that machine learning methods can map. In this study, we limit the solutions to within one orbital revolution and a semi-major axis range defined by the position geometries. To further preprocess the data, we introduce a geometry-aware normalization based on the r i-r t position pairs. With this transformation of the data, we introduce properties to the solution where it tends toward a constant value as the target radius r approaches r max (10 CDU) or the target angle θ approaches π. These properties can then be linked to structures in the trained DeepONet and KAN models. However, a comprehensive Lambert solver must accommodate a broader range of orbital configurations beyond single elliptical orbits (Russell, 2022; Negrete and Abdelkhalik, 2025). Extending this geometry-aware normalization methodology to address more complex solution spaces remains an avenue for future research.

Within the normalized formulation, the MLP provides the strongest predictive performance, while DeepONet and KAN offer more structured parameterizations of the same inverse map. DeepONet separates the dependence on transfer geometry from that on time-of-flight, allowing the solution to be analyzed in terms of branch coefficients and trunk functions. KAN, in turn, represents the mapping through learned edge functions that reveal how the solution varies with r, θ, and τ. These structured models provide complementary views of the normalized Lambert map, although neither matches the MLP’s predictive accuracy.

The mission-design case studies show that the learned surrogates can preserve the broad structure of the transfer landscape and identify promising transfer regions, which is the primary requirement in preliminary trajectory analysis. As expected, the classical iterative method provides higher local accuracy, since it is designed to converge to the required solution tolerance. The results further indicate that classical Lambert solvers remain preferable for smaller problems, whereas the neural surrogates can offer modest computational advantages in large-scale evaluations, at the expense of increased solution error. Taken together, these findings support the use of geometry-aware neural surrogates in early-stage mission design, while also indicating that improved robustness, accuracy, and extension to broader orbital configurations remain important directions for future work.

7 Conclusion

In this work, we studied the bounded zero-revolution elliptic Lambert problem and showed that a geometry-aware normalization can regularize its inverse map for neural surrogate modeling. By representing each transfer instance in canonical coordinates with bounded, geometry-dependent ranges, the proposed formulation brings diverse transfer geometries into a common form and makes the learning task more tractable. We focus on learning the inverse function in the Lagrange formulation, enabling a more structured representation of the solution space through the use of three different neural network architectures, MLP, DeepONet and KAN. The neural networks were trained to predict the semi-major axis of the transfer ellipse from the normalized variables. Among the models tested, the MLP achieved the highest predictive accuracy.

More broadly, the results suggest that the main contribution of the approach lies in the normalization itself, which enables lightweight surrogates to capture the dominant structure of the bounded elliptic Lambert solution space. We further show that this structure, when combined with appropriate network architectures, reveals a clear trade-off between inductive bias and approximation error. In Earth–Mars and multi-planetary flyby examples, the surrogates preserve the broad features of the classical solutions and identify promising transfer regions. This supports the use of geometry-aware normalization and neural surrogates in preliminary mission design and early exploration of transfer opportunities.

Generative AI statement

The author(s) declared that generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/frspt.2026.1772629/full#supplementary-material

Supplementary material Supplementary material at the source

Footnotes

Edited by: Ram Krishan Sharma, Karunya Institute of Technology and Sciences, India

Reviewed by: Yangyang Ma, Northwestern Polytechnical University, China

William Brandenburg, The University of Texas at Austin, United States

Footnotes

The target position will be passed each orbit, but the solution will ensure that the spacecraft is at the target position at t t = t i + τ.

https://github.com/KindXiaoming/pykan