The Koopman Operator

The Koopman Operator (KO) can solve for Halo orbits. The figure reports two sets of 30 different orbits calculated through a KO using basis functions of order 5, around L1 and L2. The figure includes the projection of the orbits in the three reference planes and an isometric three-dimensional view of the trajectories. The black point in the picture represents the secondary celestial body, which is Earth in the selected Sun–Earth system. As the figure shows, the Koopman solution achieves a robust representation of the orbits, the accuracy of which will be assessed later. The periodicity is given by an accurate selection of the initial conditions, which have been found through the implementation of a differential corrector in such a way that the proposed solution describes the same orbits for multiple revolutions. The KO is able to accurately represent orbits around both libration points.

The eigenfunction analysis and the study of the minimum points (minima) are important when seeking stability and periodic orbits in a chaotic application, such as the CRTBP. The KO orbits are able to maintain their periodicity due to the eigendecomposition of the dynamics, where the system is diagonalized and each eigenfunction is studied in terms of magnitude and shape. Therefore, if desired, unstable contributions from eigenfunctions connected to unstable eigenvalues can be separated. The figure reports the first 24 eigenfunctions from the third-order KO solution. The six-dimensional functions are represented with just two dimensions in the (x; vy) plane, because this information sufficiently describes each Lyapunov orbit. The eigenfunctions are centered at the L1 point, meaning that position (0, 0, 0) corresponds to L1 instead of the primary celestial body. Examining a gray-scale gradient, the figure shows the locations of minima in bright white. These curves represent equilibrium conditions, or more precisely, their approximation through a third-order KO. Therefore, when we include the initial conditions in the (x; vy) plane as white dots, their locations settle close to the minimum lines. Thus, the Koopman approximation is correctly approximating the dynamics of the system in terms of its eigenfunctions. However, the accuracy of the approximation decreases as the initial condition becomes larger; the farthest white points are not as well approximated as the white points near the origin. Moreover, the location of the minima is not unique. Many eigenfunctions show well-marked bright lines in the bottom part of the plane, for low velocities. These initial conditions are still connected to periodic orbits. The resulting pathway leaves the Lagrangian L1 point and starts orbiting around the main celestial body with a periodic behavior. Therefore, even if the eigenfunctions become less accurate the farther they are from the equilibrium point, their minima are still able to connect to different kinds of periodic orbits. Note also that the figure reports only initial conditions for orbits connected to negative x and positive vy. However, thanks to the symmetric properties of the CRTBP, the domain of the eigenfunctions can be widened to support retrograde orbits and all the remaining periodic orbits from axial and central symmetries.

A spectral study can also be performmed using KO. The figure compares the eigenvalues, computed using the Galerkin method, from the third-order and sixth-order Koopman matrices. In this particular application, the eigenvalues of the third-order solution are a subset of the sixth-order approximation. Thus, the sixth-order solution incorporates information of previous orders and enhances it with additional nonlinear functions. As such, the location of the high-order eigenvalues is at the edge of the diamond-shaped figure. By looking at the intensity of every singlemarker, it can be noted how eigenvalues closer to the origin and the imaginary axis are repeated multiple times. Indeed, it is possible to know in advance the position of the unperturbed eigenvalues, for any given order, by the simple composition of the eigenvalues that have already been calculated for lower orders. Therefore, position and repetitiveness are explained in terms of statics and combinatorics (the theory of combinations), where all the new eigenvalues are found as a linear combination of the original ones. Moreover, the linear approximation near the L1 is known analytically, aswell as its eigenvalues,which have been portrayed in the figure as circles. The linear part of the KO solution matches the linear contribution of the true dynamics. The six linear eigenvalues are perfectly represented, and the high-order approximations evolve from the latter to best represent the nonlinear terms of the system of ODEs.

The frequencies and spectral behavior of the system provided by the KO can also be used to assess the convergence of the methodology. This can be done by the study of theKoopman modes, that is, the projection of the observables in the operator.The figure presents the intensity of the Koopman modes of the system. Particularly, it shows how the eigenvalues with the highest modulus are the least influential in the solution. However, the small contributions from the modes connected to large eigenvalues are crucial to achieving an accurate approximation of the dynamics, and they make the difference in accuracy when comparing the order 3 solution with the one of order 6. Therefore, many terms other than the linear ones are nonzero, with influence that decreases as the order of the eigenfunctions increases. This also implies that the Koopman modes decay with larger frequency values, leading to a possible truncation of the expansion for eigenfunctions with particularly high frequencies.

The figure shows, with different line styles, the accuracy of the approximation on each orbit for different orders of the KO. As expected, as the halo orbits become larger, the position error increases consequently. Therefore, for any given order, the position error convergence analysis describes a set of curves. In that regard, order 3, represented by a continuous line, is the least-accurate solution, while order 6 is the most accurate. However, it is worth noticing how, for each orbit, the gain in accuracy obtained by increasing the KO order behaves differently depending on the initial condition. The sixth-order KO solution is two orders of magnitude more accurate than the third order for orbits close to L1. However, as the initial condition of the halo orbit moves away from the equilibrium point, the gain in accuracy is reduced, settling around a single order of magnitude. Indeed, for very large initial conditions, increasing the order of the basis functions to over a given value stops being beneficial because the gain in performance does not repay the higher computational burden. This aspect is connected with the Richardson approximation applied for the polynomial representation of the dynamics, which fails to provide an accurate approximation of the system for initial conditions far from the equilibrium point. Consequently, it can be noted that for each set of curves, at different orders, there is a well-marked separation in the bottom part of the graph, while curves from different orders start to get closer as the orbit becomes larger. This aspect has been highlighted in the right figure, where the mean position error  for each orbit is compared with the size of the halo orbit and the maximum order of the basis functions used to represent the KO. Thus, orbits with the furthest initial conditions from the L1 point generate mean position errors that range in value by up to one order of magnitude. However, as we get closer to L1, the mean position error for each KO solution order considered decreases with a different slope, making the accuracy gain more substantial as we increase the order of basis functions used. Additionally, it can also be observed that there is a region very close to L1 where the error seems to reach convergence and stabilize for six order basis functions. Thus, for halo orbits close to L1, an increase of the KO order leads to a substantial improvement in accuracy and in the performance of the approximation of the system.

A new particle filter based on the KO is proposed, where each integration is substituted by the mere evaluations of functions using the KO approximation of the solution for the dynamics. In particular, the Koopman operator particle filter (KOPF) evaluates the propagated particles, each time step, by considering their different initial conditions during the evaluation of the eigenfunctions. After the dynamics have been analyzed and linearized in the set of basis functions, the state of each particle, for any given time, is a function of its own initial conditions. Each particle is propagated through the evaluation of different initial conditions into the eigenfunctions. Therefore, all the particles are integrated in the KO framework, and the state probability density function is propagated rapidly, making the KOPF suitable for onboard applications.

The Koopman theory is also integrated into the control segment of the station-keeping application. The KO solution is therefore merged inside a model predictive controller, where the control feedback input is optimized by predicting the future behavior of the state of the system given the current state. The basic NMPC algorithm can be summarized in three major steps. At first, the state of the system is measured. This step is achieved through the state estimation performed by the KOPF. Then, an optimal control problem is solved, with the goal of minimizing some well-defined cost function with respect to the control input and subject to either control or state constraints. Lastly, the NMPC feedback value is selected and used as a control value for the next sampling period. The KO is used inside the NMPC for the prediction horizon length, where the state of the system is propagated forward in time to assess the validity of the proposed control input value. Therefore, the newly proposed KO solution of the CRTBP is embedded both in the estimator and in the controller part of the closed-loop station-keeping algorithm.

An overlook of the system architecture of the KO uncertainty propagation technique and of the KOF is proposed in the figure. The KO decomposition of the system, in terms of eigenfunctions and observable matrices, is performed offline: the dynamics are fed into the KO framework, which outputs the state transition polynomial map (STPM) that propagates the state for a specific time length. That is, for onboard navigation, it is sufficient to just upload the STPM on the onboard computer because the evaluation of the inner products and the evaluation of the observable matrices need to be computed only once. After storing the polynomials, the filtering algorithm works autonomously, predicting the state uncertainties and processing measurements. First, the STPM is shifted to be centered at the current estimate, and the Isserlis moments are evaluated given the current covariance. The state PDF is then propagated via its central moments through the evaluation of the shifted STPM, where the Isserlis moments are substituted to their monomials counterparts. The KO uncertainty propagation technique can approximate high-order central moments up to any arbitrary order for analysis purposes or to fit a polynomial update. Figure 1 shows the simple case of a linear measurement update, where just the predicted state mean and covariance are used to evaluate their corrected counterparts. Thus, the KOF prediction step is completed by the evaluation of the predicted measurement mean and covariance in the KO framework. Indeed, the measurement equations have been analyzed offline to evaluate the measurement KO observable matrix, which defines a new polynomial map that directly connects the state variable to the measurement. After performing the polynomial shift such that the shifted map connects state deviation vectors to the measurement outcomes, the measurement moments are likewise predicted for the state. Lastly, when the true measurement becomes available from the sensors, the update is computed and the filter can start the process over with a new prediction step. The system architecture shows how the KO methodology derived in previous works has been exploited and modified to work with stochastic variables and to perform recursive estimation. The STPM shifting and evaluation are critical for the correct prediction of the moments. Moreover, the measurement can be fully represented in the KOframework as well such that it is possible to create a measurement polynomial map that undergoes the same calculation of expectations.

A Monte Carlo analysis has been performed to assess the consistency of the KOF. Figure 2 shows the state estimation for 3000 runs, with random true initial condition drawn according to the initial Gaussian distribution. The figure shows the errors of the position (first row) and velocity (second row) of each run in gray. The mean of the errors is reported as a black line. This mean line settles on the zero value rapidly after every restart of the oscillator, which is in accordance with the unbiased nature of the Koopman operator filter. The convergence is assessed by the overall decrease of the uncertainties of the state of the system. The initial covariance is particularly large, and the figure zooms in to show the filter behavior at steady state, where the oscillations are marked by the sinusoidal pattern of the covariance lines. Indeed, the blue lines represent the effective standard deviation of the error as 3σ, which are calculated directly from the Monte Carlo runs at each time step; whereas the cyan lines are the estimated uncertainties from the updated covariance matrix of the filter, which are plotted again as 3σ. A consistent filter has a correct prediction of its own spread of the error, which is assessed in the figure by the overlapping of the estimated covariance (cyan) over the effective covariance (blue). The effective standard deviation indicates the effective level of accuracy of the filter, and it is evaluated as a mean among the runs. This parameter directly evaluates the statistics of the errors from the results of the Monte Carlo analysis, and it is used for validation purposes. On the other hand, the predicted standard deviation indicates the level of accuracy the filter is estimating it is reaching and it is evaluated directly from the updated covariance matrix of the filter. Thus, a consistent filter is able to correctly estimate its own accuracy level, and there is a match between the effective and predicted covariance. For the Duffing oscillator, the KOF has been proven to be a coherent and consistent filter, without any estimation bias; and it is robust against brute restarts of the system.

The newly developed technique has been applied to the CRTBP in the Earth–moon system, and the results have been validated through a Monte Carlo analysis,where the predicted  moments of the state have been compared to the effective moments evaluated from all of the Monte Carlo runs.  The figure shows the time evolution of the uncertainties represented with the state raw moments. Indeed, the continuous lines represents the true evolution of the moments of the system: the mean is in black; the covariance, as the 3σ standard deviation, is in blue; the skewness is in cyan; and the kurtosis is inmagenta. The correct prediction of the  mean and uncertainties of the systemis asserted by the overlapping of the points over the continuous lines.The accuracy of the estimation of the moments decreases as the integrating time increases, meaning that, for extremely large propagation times, the prediction of the state PDF is unfeasible. The CRTBP is strongly diverging: thus, high-order central moments after long time steps require an elevated amount of runs in order to have correct indicating values. Moreover, the figure shows the high level of accuracy in the first half of the simulation for times shorter than one year, which is zoomed in on in the figure. Stationkeeping applications, in filtering, usually receive measurements in the section where the KO shows an extremely accurate prediction ofmoments. The KO technique can predict central moments up to any arbitrary high order. However, due to the fact that odd moments can be negative, the accuracy of the prediction of even moments is higher than their odd counterparts. Indeed, some values of the skewness, in the figure, are slightly off by the end of the simulation, where the state PDF largely expands without control. In the remaining sectors, the prediction of odd moments outperforms any other uncertainty prediction technique based on the Gaussian approximation of the state distribution (such as UT) because they provide a null value for any odd moment.

The central moments estimation can be appreciated by reporting the uncertainties as around their mean. Indeed, although the previous figure shows the correct prediction of the family of orbits in terms of its overall possible trajectories, this figure focuses on the spread around the mean. Therefore, this figure displays similar results as the previous one, but with a new prospective centered at the current mean. It shows the chaotic nature of the problem because the uncertainties, in terms of their σ value, increase exponentially; and the spread of all the possible resolutions of the state creates a conical shape. The KO prediction is able to keep track of the central moments of the system, especially during the first half of the simulation. Both the state covariance and kurtosis expand without any boundaries. The star points show the correct estimation performed by the KO technique. The new methodology is keeping an accurate prediction of the state distribution, even for a long period of time, and even if slightly reducing the accuracy. The figure has a logarithmic scale for the values of the moments, highlighting the diverging nature of the CRTBP dynamics. Thus, it is worth noticing the exponential increase of the spread of the state PDF, with standard deviations that are two orders of magnitude bigger after merely three-quarters of a full revolution.

The correct prediction of the uncertainties has been compared to other common techniques that propagate central moments forward in time. First, the Gaussian mixture model has been used to propagate the state distribution by splitting the PDF in 2n + 1 Gaussians. Then, state transition tensors and differential algebra have been adopted to provide a high-accuracy level comparison of the state propagated standard deviation. These two techniques have been proven to obtain similar results and the same level of accuracy when applied up to a high truncation order. The figure reports the percentage relative error level with respect to the Monte Carlo standard deviation values of the position ζp and velocity ζv covariance prediction. The GMM approximation has difficulties in estimating the spread of the distribution and the line relative to its performance (in yellow) is far above the STT and DA lines (in red) and the KO line (in blue). By zooming in on the lower error magnitude, it can be noted that the high-order approximation with STTs and DA has a smaller relative error when compared to the KO. This aspect is due to the Richardson approximation of the dynamics of the CRTBP, which becomes less accurate for halo orbits as the time increases. The KO technique has evaluated its state transition polynomial map on the Richardson approximation of the Hamiltonian; for this specific application, it suffers the same limitations.

The figure shows, for each component of the state, the Monte Carlo analysis performed with the KOF on a selected Lyapunov orbit. Indeed, each gray line represents the errors connected to one single run of the Monte Carlo, for each state component. The means of the errors are portrayed with black lines. The mean (black) lines settle on the zero value, proving that the KOF is an unbiased filter, as expected from the linearized minimum mean square error principle on which the update of the filter is based. The Monte Carlo analysis provides also information about the spread of the state error and the level of uncertainties. The figure has multiple pairs of dashed and continuous blue lines. The dashed blue lines represent the levels of standard deviations, which are obtained from the Monte Carlo runs: they represent the effective level of accuracy of the filter. These effective lines are calculated, for each time step, by considering the values of the state errors among all the simulations; and they represent the actual behavior of the filter. On the contrary, the continuous blue lines report the estimated level of the uncertainties, which is predicted by the filter. Thus, for each time step, these lines are calculated by extracting the error standard deviation directly from the updated covariance of the filter. By looking at the figure, we can assess the consistency of the KOF because the effective and predicted lines overlap, meaning that the filter is able to correctly predict its own uncertainty. The KOF shows convergence and consistency, where steady-state accuracy levels are reached rapidly after a few updates. The large initial uncertainties are rapidly reduced around the most current estimate, and the updated covariance of the filter correctly represents the spread of the state error PDF.

The figure reports the accuracy levels, in terms of error standard deviations, of four filters, displaying estimated values with continuous lines and the actual behavior with dashed lines. The EKF is shown in red. The effective EKF lines, for both the position and velocity, are out of scale when compared to the predicted counterparts, which are overlapping with the predicted IKF lines. Indeed, the green lines are connected to the IKF. The EKF and the IKF share the same updated covariance matrix while their estimate changes. Thus, the overlapping between the IKF and the EKF predicted lines matches with what is expected from theory: the two estimators approximate the posterior distribution of the state as Gaussian and with the same level of uncertainty (same standard deviations) but different means. Indeed, although the EKF is a minimum mean square error (MMSE) filter, the IKF is a filter based on the MAP principle, which outputs as its estimate the most likely state of the posterior distribution. However, the IKF and EKF share the same prediction step, based on the linearization of the dynamics, which is not sufficient to achieve an accurate approximation of the state prior distribution. Therefore, both filters are inconsistent; and their dashed lines settle orders of magnitude above the continuous ones. The EKF and IKF believe that they are achieving higher accuracy than their actual results, and they are overconfident in their performance. The difference between the green dashed lines and the red dashed lines is connected to the different update steps performed by the two filtering techniques: a MMSE one for the EKF, and a MAP one for the IKF. On the contrary, in blue, the UKF applies the unscented transformation in its prediction step to obtain a more accurate state prior distribution. Thus, the effective UKF line settles below the EKF and the IKF ones. However, the matching with the predicted uncertainty is missing and the filter is inconsistent as well. The UKF estimates a performance similar to the KOF (in black), which is shown by the overlapping between the predicted KOF and the predicted UKF lines. Therefore, the KOF is the only filter that converges with consistency. The KOF is able to correctly predict its own uncertainties, and its estimate is the most accurate among the other filtering techniques. The effective KOF lines are the lowest effective lines for both position and velocity, and they overlap their predicted counterparts along the whole time length of the simulation.

The figurs shows the Duffing Oscillator system brought to the origin for different amounts of times. The uncontrolled system has the mass oscillating between −1 and 1 indefinitely, displayed in black. The energy-optimal approach has the state of the system rotating in a spiral around the equilibrium point, at a rate that gets slower, the larger the final time. Thus, the required magnitude of control effort decreases as the system is allowed more time to reach stasis.

A second analysis is performed with a fixed required time of  2 seconds, where the initial condition of the state is picked following a grid, and the system is lead to stasis. The figure shows the path to equilibrium from the 12 different initial conditions, and the relative uncontrolled behavior in the background. Therefore, control has been achieved
given different states and time constraints.

The controlled pathway of the chaser spaecraft to meet the target spacecraft can be appreciated in the figure, where it can be noted that the spacecraft slows down and changes its direction to approach the origin with smaller oscillations each time. The velocity profile, the right plot of the figure, shows how the chaser decreases its velocity in a spiral, reaching stasis at the desired final time. The rendezvous is achieved and the control technique has been proven to work successfully.

The figure shows an application that aims to achieve rendezvous when the chaser is 10 km ahead of the target, given the constraint of performing the maneuver in 12 hours. The figure shows the resulting pathways of the chaser subject to different dynamics. The uncontrolled linear propagation results to a constant integration, which is a single point in the figure. Therefore, the uncontrolled linear solution, in green, is just a point because it is static. On the contrary, the uncontrolled nonlinear dynamics, in black, has been evaluated considering the first six potential polynomials. This line shows the actual behavior related to the uncontrolled relative motion of the spacecraft, with the presence of the translational drift. The controlled pathways are reported in the figure with red and blue lines for the linear and the high-order dynamics, respectively. Thus, the correct solution of the controlled chaser can be appreciated by analyzing the blue curves. As the chaser gains velocity by increasing its x components, it is able to rapidly approach the target in the y direction. After passing the half-point, the control starts to slow the chaser down, such that it reaches rendezvous with null velocity. The controlled linear pathway, in red, is reported for comparison purposes. The red curve shows the behavior of the chaser that is propagated under linearized dynamics, which fails to achieve rendezvous when controlled by the accurate nonlinear control law. Thus, the consideration of nonlinear terms is crucial for the high accuracy and precision evaluation of the costates of the system, and the newly developed energy-optimal control technique is able to evaluate such control.

The rendezvous application has been repeated with the chaser’s initial condition at different relative positions, in order to assess the capability of the KO technique to evaluate the optimal control to approach the target from every direction. As such, the figure reports the optimal rendezvous pathways of the chaser starting from a distance of 2 km from the target (dotted circle), with initial condition separated every π∕4. The figure reports solutions for rendezvous in exactly 4 h. The values of the initial relative velocity are evaluated considering the difference in the orbit’s tangential velocity, where the reference orbit radius a is increased by a factor x. This aspect best describes the relative initial condition between chaser and target orbiting in different orbits. Therefore, the figure shows that the KO energy-optimal inverse control technique is able to calculate the exact costates for the chaser approaching from every direction. These solutions have been calculated using a high-order Lagrangian that considers the first six polynomials in the expansion of the potential energy.

Solving Lambert’s Problem using multiple revolutions involves finding transfer orbits that complete one, or more, revolutions around the gravitational body before reaching
the desired final position vector with a given time of flight. A KO polynomial map is created for each number of revolutions without the need to recalculate the KO matrix, a benefit from using the KO theory. The possible solutions to the multiple revolutions Lambert’s problem are shown in the figure for a variety of transfer times and semi-major axes. This figure presents the comprehensive solution space, encapsulating both direct transfer orbits and transfer orbits that complete multiple revolutions. It offers an entire view of the potential transfers a satellite can take between the specified position vectors. There are seven possible orbital trajectories for transfer times greater than about 4800 s: two possible multiple revolutions transfer orbits for 𝑁 = 1, 2 and 3 and one possible solution for a single revolution orbit transfer. Such times are too long for another single revolution orbit transfer solution as this other possible solution approaches the time corresponding to a parabolic transfer orbit for increasing semi-major axis. The 0, 1, 2, and 3 revolutions solutions are plotted in the left figure. These correspond in color to the multiple revolutions solutions plotted in the right figure. Some of these multiple revolution solutions are not feasible due to having small semi-major axes that cause their orbital trajectories to pass through the Earth for the time constraint considered. Using the KO methodology allows for the calculation of multiple revolutions solutions by increasing the number 𝑁 without having to change the time of flight. This entails benefits because the zonal representation uses true anomaly as its independent variable.

By working with true anomaly as the independent variable, the Koopman Operator State Transition Polynomial Map (KOSTPM) finds the curve, expressed as interpolated polynomials, in the velocity space of orbits that are a candidate for Lambert’s problem. The curve in the  velocity domain and the relative family of orbits are reported in blue, respectively. Every solution is propagated using the same Koopman matrix that connects the position vectors in the TPBVP, as the change in the anomaly of the transfer orbit is the same regardless of the transfer time. That is, looking at the family of transfer orbits, the left figure, they all have a different semi-major axis, a different time of flight, and a different velocity vector, but they all cover the same true anomaly span, Δ𝜃. Therefore, for a specific point in the spacecraft orbit, the family of velocities that create orbits connecting the two points is evaluated in the KO framework. This curve, called Ξ(𝜃), is represented in the three-dimensional velocity space in the right figure. The analytical function of the curve suits the proposed methodology for the application of various optimization techniques and analysis. For example, since each orbit starts from the same position vector, i.e., constant range 𝑟, the orbit transfer with the minimum energy is the one with the minimum velocity. The orbit with the minimum energy is connected to the hyperbola’s shortest velocity vector. The figure highlights the point of Ξ(𝜃) closest to the origin, which corresponds to the minimum velocity norm and, therefore, to the minimum energy transfer. The transfer is highlighted in red. At the same time, it identifies the smallest orbit with the lowest semi-major axis since it is the one with the lowest specific energy value.

The family of possible transfer orbits has been parameterized for each point of curve Ξ(𝜃). This representation leads to other kinds of optimizations, such as minimum fuel transfer, with the assumption of impulsive maneuvers. Thus, assume that the spacecraft is moving between two separate orbits rather than two points in space. Knowing the velocity vector at the giving point and the family of transfer orbit Ξ(𝜃), the minimum fuel solution corresponds to the transfer orbit that requires the minimum total impulse to let the spacecraft change from initial to transfer orbit and from transfer to final orbit. The left figure reports in red the minimum fuel solution, while the right one shows the impulse vectors in the velocity space as the difference between the velocity hyperbolas at the initial and final conditions. The minimum condition is obtained with a global search on the sum of two impulses to identify the transfer orbit that requires the smallest overall change in velocity. Knowing the velocity hyperbola at the initial and final position, thanks to the KOSTPM leads to an easy scalar optimization problem among vectors’ magnitudes.

Time optimization can be performed under the assumption that a maximum amount of fuel is available. Assuming a maximum total impulse of Δ𝑉𝑀𝐴𝑋 = 3 km/s, the minimum time of flight optimization is translated to determining the transfer orbit with a total impulse equal to the constraint, as that corresponds to the fastest transfer orbit. That is, all the fuel is used to achieve the maximum change in velocity to set the spacecraft to the fastest route. The figure  shows the maximum impulse constraint on the total impulse curve as a red horizontal line. The constraint crosses the function in two separate points, Δ𝑉1 and Δ𝑉2. One point corresponds to the minimum time-of-flight solution, while the other is its counterpart, i.e., the orbit with the same change in velocity but selected on the opposite side with respect to the minimum impulse orbit. The solution pathway is reported as well, where the light blue orbit is the optimal time orbit, while the red orbit is the minimum impulse one from the previous analysis. The magenta orbit is the counterpart.

It is possible to calculate the minimum Δ𝑣 trajectory given any initial and final position vectors on two generic elliptic orbits. Hence, it is possible to produce “pork-chop” plots that identify the best departure and arrival orbital positions on these orbits. The minimum Δ𝑣 direct transfer orbit connecting each possible initial true anomaly and each possible final true anomaly along the elliptic orbits is calculated. Values are interpolated so that the minimum Δ𝑣 is calculated every 1 deg in true anomaly. The resulting matrix is displayed as a contour map. A generic pattern can be seen in the figure where diagonals across the matrix have similar Δ𝑣 values. The difference between the true anomaly at arrival and the  true anomaly at departure is constant along these lines. The Δ𝑣 required to transfer between these true anomalies is a significant factor in determining the overall Δ𝑣 requirements for the trajectory. The constant true anomaly difference along each diagonal results in similar minimum Δ𝑣 values for each trajectory on that diagonal. A 3-D “pork-chop” plot of the matrix plot is shown as well. The 3-D pork-chop plot and the contour plot can be used to identify the best departure and arrival orbital positions that correspond to the overall minimum Δ𝑣 transfer.

The distribution of eigenvalues of the system for the solution that uses order 9 in the basis functions is shown in the figure. Here the maximum values of the imaginary part of the KO eigenvalues is 9, respecting the duplicity of the eigenvalues creation. The overall shape and location of the eigenvalues is kept ifor different KO orders. Hence, the overall spectral decomposition of the dynamics remains, however, it becomes more accurate as the order is increased which increases the number of eigenvalues that are farther from the origin. This eigenvalue distribution is representative of the whole J2 dynamical system for any initial conditions that are considered for solving Lambert’s problem since the Koopman Matrix is independent of these conditions. To better understand how the solution to Lambert’s problem in the perturbed system depends on each Koopman mode, the Koopman mode magnitudes are plotted against the absolute values of the Koopman eigenvalues. The magnitudes are plotted for orders 7, 8, and 9 in the KO basis functions. The number of modes increases with an increase in order, however the magnitudes of these modes decrease in value. Hence, increasing the order of the KO basis functions still leads to an increase in the accuracy of the solution since influential modes are being added. However, as the order is increased, the benefits to improving the accuracy of the solution decreases. When using higher-order basis functions, there exists a trade-off between the increase in computational complexity and the improvement in the accuracy of the solution.