Giant constellation distributed orbit determination method implemented on domestic super computer and solver system thereof

Abstract

The invention discloses a giant constellation distributed orbit determination solver and method based on a domestic super computer. The solver adopts a'base layer-algorithm layer-application layer-interface layer 'modular architecture, and a C-API is externally provided by an external' C 'handle, and can be directly called by Fortan / Python; a kernel realizes integration-correction integration through RKF7 (8) adaptive integration, Armijo non-precision line search and mixed precision PCG, and in combination with a CSR-D sparse format and topological perception domain decomposition, the memory complexity is reduced from O (N) to O (N) under the scale of ten thousand kernels, and the communication synchronization frequency is reduced by more than or equal to 70%. The calculation process comprises the steps of step length correction driven by local truncation errors, one-time completion of Armijo backtracking global stipulation, float-double residual error triggering type precision upgrading and Allan variance sliding window dynamic PPD early warning. The method supports check point continuous calculation and cross-language integration, and can be widely applied to operation control and space traffic management of giant constellations such as Beidou and Kiaji.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of spacecraft precise orbit determination, and in particular to a distributed orbit determination method for a giant constellation implemented on a domestic supercomputer and a solver system thereof. BACKGROUND

[0002] With the rapid deployment of super-large low-earth orbit constellations such as Starlink and OneWeb, the number of satellites has increased from hundreds to thousands or even tens of thousands, and the traditional single-machine or GPU cluster orbit determination scheme has encountered bottlenecks in terms of computing scale, memory capacity, and cross-node communication. Existing technologies (such as CN117776601A) propose a distributed least squares orbit determination framework based on MPI+OpenMP, but still have the following shortcomings:

[0003] Fixed step numerical integration and precise line search are used, resulting in frequent global reduction communication in the multi-core scenario, and the parallel efficiency decreases sharply with the increase in the number of cores;

[0004] The Jacobian matrix is stored in a dense block, and the memory occupancy is proportional to the square of the number of satellites, so the single-node memory is quickly exhausted;

[0005] Only C++ template interfaces are provided, and the legacy Fortran / IDL task system needs to be additionally encapsulated, and the cross-language call delay is high;

[0006] There is a lack of checkpointing and error positioning mechanism for long-time tasks, and when the supercomputer job fails, it needs to be calculated from the beginning, wasting a lot of machine time.

[0007] Therefore, there is an urgent need for a new giant constellation orbit determination method with low communication volume, controllable memory, easy integration across languages, and adaptation to domestic supercomputer architecture. SUMMARY

[0008] The present application aims to provide a giant constellation distributed orbit determination method implemented on a domestic supercomputer and a solver system thereof, thereby solving the aforementioned problems in the prior art.

[0009] To achieve the above-mentioned purpose, the technical solution adopted by the present application is as follows:

[0010] A constellation orbit determination solver based on a domestic supercomputer, the solver is implemented with a standard C++ kernel, and C language style function interfaces are exported externally through extern "C", and the interfaces manage the whole life cycle of the solver in handle mode;

[0011] The solver is divided into four modules: a basic layer, an algorithm layer, an application layer, and an interface and management layer, wherein:

[0012] The basic layer includes:

[0013] a) Distributed sparse matrix storage with Eigen::SparseMatrix<…> and support for CSR, CSC, COO formats;

[0014] b) Memory pool management operator based on memory pool pre-allocation + object reuse strategy, operator calls vector.reserve() to pre-allocate capacity before orbit sequence construction and reuses stack buffer during interpolation process;

[0015] c) Earth non-spherical gravity perturbation model, at least containing J2, J3, J4 zonal terms, and reserving spherical harmonic coefficient expansion interface;

[0016] d) Space-time reference operator, realizing geodetic system and geodetic inertial system coordinate conversion and absolute calendar time and relative second exchange;

[0017] Algorithm layer includes:

[0018] a) Orbit calculation operator, using RKF7(8) adaptive step-size integrator, integration process is completed completely in geodetic inertial system, step-size control through 7th and 8th order solution difference value estimation local error;

[0019] b) Linear system solving operator, for equation AT·A·Δx=–AT·r, automatically switching according to non-zero element proportion threshold: when sparsity ≥ 95%, using conjugate gradient method, otherwise using LU decomposition;

[0020] c) Differential correction operator, integrating non-precise line search based on Armijo condition, step-size α initial value taking 1.0 and decreasing in backtracking way until satisfying sufficient descent condition;

[0021] Application layer includes:

[0022] a) Orbit determination operator, transforming giant constellation orbit determination problem into “prediction-comparison-correction” closed loop iteration, iteration convergence criterion being consecutive two times residual RMS improvement amount less than set threshold;

[0023] b) Precision evaluation operator, calculating parameter covariance matrix based on normal equation inverse matrix and propagating to obtain position precision factor PPD and velocity precision factor VPD;

[0024] c) Visualization operator, using command line single line return cover technology to output progress bar and Statistics online statistical summary in real time, Statistics class outputting mean, standard deviation with O(1) space complexity;

[0025] Interface and management layer includes:

[0026] a) The programming interface is encapsulated, and all core functions are provided in the form of C functions prefixed with orbit_solver_. Error codes are returned using the OrbitSolverStatus enumeration.

[0027] b) Management interface adopts a hybrid parallel architecture of MPI+OpenMP. MPI distributes the satellite constellation to different computing nodes according to the orbital plane. OpenMP parallelizes the propagation of single-satellite orbits within the node and realizes parallel file reading and writing through MPI-IO.

[0028] c) Exception handling interface, providing get_last_error() to obtain the error code, location of occurrence and possible cause, supports checkpoint setting and recovery mechanism, and enables the continuation of long-term task interruption;

[0029] The solver uses a domain decomposition method to distribute the sparse Jacobian matrix formed by the joint orbit determination of tens of thousands of satellites to domestic supercomputing nodes with tens of thousands of cores, thereby reducing the overall memory usage and computational complexity by an order of magnitude. It can also converge the initial orbital error of hundreds of meters to the centimeter level within 10-15 iterations.

[0030] Furthermore, the distributed sparse matrix operation framework achieves communication hiding and storage optimization in the following ways:

[0031] The Jacobian matrix is ​​divided into blocks and rows according to the orbital plane-satellite number to form a submatrix A_i that corresponds one-to-one with the supercomputing node topology;

[0032] A double buffering mechanism is adopted. While performing sparse matrix-vector multiplication on the local submatrix A_i, the neighbor boundary vector Δx_{i±1} is exchanged through the MPI non-blocking interface Iallreduce, so as to achieve overlap of computation and communication.

[0033] Use MPI_Dist_graph_create_adjacent to build an adjacent communication graph of the orbital plane, fix the communication degree at 2, and reduce network congestion;

[0034] Based on the orbital dynamics Jacobi sparse structure, a dedicated CSR-D format is used for storage: row pointers only record satellite rows containing non-zero elements, and column indexes are stored in segments according to perturbation source type to reduce index memory and improve cache hit rate;

[0035] The analytically estimated block diagonal approximate inverse matrix M^{-1}_i is used as the preconditioner for domain decomposition to replace the dense Schur complement, thereby reducing the memory complexity of a single node to O(N).

[0036] Furthermore, when employing the conjugate gradient method, the linear system solver uses incomplete intra-node LU decomposition and inter-node domain decomposition preconditioners, and supports mixed-precision iteration, dynamically improving accuracy in the later stages of iteration to ensure numerical stability.

[0037] Furthermore, the visualization operator dynamically refreshes the progress percentage within a single line on the terminal using the ProgressBar class provided by progress_bar.hpp, and combines it with the Statistics class to provide real-time statistical output of satellite drag coefficient, orbital altitude, and residual RMS.

[0038] Furthermore, the management interface dynamically adjusts the load of the MPI process based on the satellite observation density. When the observation data of a certain orbital plane suddenly increases, the main process re-divides the satellite list and issues a new domain decomposition scheme to achieve dynamic load balancing.

[0039] Furthermore, the exception handling interface automatically triggers a checkpoint write when it detects a std::bad_alloc or matrix condition number exceeding the limit exception, and resumes iteration from the most recent checkpoint through the orbit_solver_restore() function when the task is resubmitted.

[0040] In another embodiment, a method for determining the orbit of a giant constellation using a decoder includes:

[0041] S1. Read the initial orbital elements, station coordinates, and inter-satellite / satellite-to-ground observation data;

[0042] S2. The basic layer operator converts the data into a ground inertial system state vector and constructs a sparse Jacobian matrix;

[0043] S3. The algorithm layer operator performs RKF7(8) integration, Armijo differential correction and PCG solution until the convergence criterion is met;

[0044] S4. Application layer operators output centimeter-level orbital data, PPD / VPD accuracy indicators, and visual statistical reports;

[0045] S5. The interface layer returns track results and error codes via C functions, and supports external Python / Fortran calls.

[0046] Furthermore, the specific operational steps in step S3 are as follows:

[0047] a) Adaptive step size control: The local truncation error Δ = ||y is calculated using the RKF7(8) integrator. n +1(7)-y n+1(8‖∞,When Δ>ε_max, the step size is reduced by h_new=h·min(5, 0.8·(ε_max / Δ)^(1 / 8));

[0048] b) Armijo Inaccurate Line Search: Along the Correction Direction d k =-Δx backtracks, iteratively updating α←β·α (β=0.5) until Φ(x) is satisfied. k +αd k )≤Φ(x k )+c1α∇Φ k ᵀd k This reduces the number of global synchronized Allreduce operations by 70%.

[0049] c) Mixed-precision PCG solution: The first 80% of iterations use single-precision float, and the last 20% of iterations switch to double. Precision upgrade is automatically triggered by the relative change in residuals, and the final residual is <10⁻. 4 m and the calculation time is reduced by ≥15%.

[0050] In another embodiment, a computer-readable storage medium stores a computer program that, when executed by a processor, implements the steps of the above-described method.

[0051] In another embodiment, an electronic device includes a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the program to implement the steps of the method described above.

[0052] The beneficial effects of this invention are:

[0053] Communication hiding: By using double buffering and non-blocking communication, computation is overlapped with global reduction, significantly reducing synchronization wait.

[0054] Memory optimization: CSR-D sparse format stores only non-zero blocks, and domain decomposition preconditioners replace dense Schur complement, thus linearizing memory complexity.

[0055] Precision closed-loop: Adaptive step size control and inaccurate line search ensure rapid convergence, hybrid precision iteration balances efficiency and numerical stability, and Allan variance sliding window enables dynamic precision early warning.

[0056] Cross-language integration: The C-API handle pattern allows upper-level systems to call it without recompiling, reducing project migration costs.

[0057] Fault-tolerant continuation: The checkpoint-error code-recovery three-in-one design supports rapid resumption of operation after job interruption, improving the utilization rate of supercomputing resources.

[0058] This invention provides an efficient, robust, and easy-to-use core tool for real-time precise orbit determination of giant constellations, and has good prospects for engineering promotion. Attached Figure Description

[0059] Figure 1 This is a schematic diagram of the solver module of the present invention;

[0060] Figure 2 This is a flowchart of the solve_CG function of the present invention;

[0061] Figure 3 This is a flowchart of the differential correction process of the present invention;

[0062] Figure 4 This is a flowchart of the initial orbital data reading process of the present invention;

[0063] Figure 5 This is a flowchart of the observation data reading process of this invention;

[0064] Figure 6 This is a flowchart of the process for reading ephemeris files according to the present invention;

[0065] Figure 7 This is a flowchart of the linear system solution operator of the present invention;

[0066] Figure 8 This is a flowchart of the differential correction operator of the present invention;

[0067] Figure 9 This is a flowchart of the trajectory determination calculation of the present invention;

[0068] Figure 10 This is a flowchart of the core steps for accuracy assessment in this invention;

[0069] Figure 11 This is a flowchart of the RKF7(8) adaptive step size of the present invention;

[0070] Figure 12 This is a comparison diagram of Armijo line search communication synchronization according to the present invention;

[0071] Figure 13 This is a diagram of the mixed-precision PCG switching strategy of the present invention. Detailed Implementation

[0072] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to the accompanying drawings. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the invention.

[0073] Reference Figures 1 to 13The example shown is a constellation orbit determination solver based on a domestically produced supercomputer. The solver is implemented with a standard C++ kernel and exports a C-style function interface through extern "C". The interface manages the entire lifecycle of the solver using a handle mode.

[0074] The solver is divided into four modules: the basic layer, the algorithm layer, the application layer, and the interface and management layer.

[0075] The base layer includes:

[0076] a) It adopts Eigen::SparseMatrix<...> and supports distributed sparse matrix storage in CSR, CSC, and COO formats;

[0077] b) A memory pool management operator based on a memory pool pre-allocation + object reuse strategy. The operator calls vector.reserve() to pre-allocate capacity before constructing the track sequence and reuses the stack buffer during the interpolation process.

[0078] c) A non-spherical gravitational perturbation model of the Earth, which includes at least J2, J3, and J4 band harmonic terms, and reserves an interface for extending spherical harmonic coefficients;

[0079] d) Spatiotemporal reference operator to realize coordinate transformation between Earth-fixed system and Earth-inertial system and conversion between absolute calendar time and relative second;

[0080] The algorithm layer includes:

[0081] a) The orbit calculation operator uses the RKF7(8) adaptive step-size integrator. The integration process is completed entirely in the Earth inertial system. The step-size control estimates the local error by the difference between the 7th and 8th order solutions.

[0082] b) For linear system solution operators, for the normal equation AT·A·Δx=–AT·r, the operation is automatically switched based on the non-zero element proportion threshold: when the sparsity is ≥95%, the conjugate gradient method is used; otherwise, LU decomposition is used.

[0083] c) Differential correction operator, integrating inexact line search based on Armijo conditions, with an initial step size α of 1.0 and decreasing backtracking until the sufficient descent condition is met;

[0084] The application layer includes:

[0085] a) The orbit determination operator transforms the giant constellation orbit determination problem into a closed-loop iteration of "prediction-comparison-correction". The convergence criterion for the iteration is that the improvement of the residual RMS is less than a set threshold for two consecutive iterations.

[0086] b) Accuracy assessment operator: The parameter covariance matrix is ​​calculated based on the inverse matrix of the normal equation and propagated to obtain the position accuracy factor PPD and velocity accuracy factor VPD.

[0087] c) Visual operators use command line single-line carriage return overlay technology to output progress bars and online statistical summaries of Statistics in real time. The Statistics class outputs the mean and standard deviation with O(1) space complexity.

[0088] The interface and management layer include:

[0089] a) The programming interface is encapsulated, and all core functions are provided in the form of C functions prefixed with orbit_solver_. Error codes are returned using the OrbitSolverStatus enumeration.

[0090] b) Management interface adopts a hybrid parallel architecture of MPI+OpenMP. MPI distributes the satellite constellation to different computing nodes according to the orbital plane. OpenMP parallelizes the propagation of single-satellite orbits within the node and realizes parallel file reading and writing through MPI-IO.

[0091] c) Exception handling interface, providing get_last_error() to obtain the error code, location of occurrence and possible cause, supports checkpoint setting and recovery mechanism, and enables the continuation of long-term task interruption;

[0092] The solver uses a domain decomposition method to distribute the sparse Jacobian matrix formed by the joint orbit determination of tens of thousands of satellites to domestic supercomputing nodes with tens of thousands of cores, thereby reducing the overall memory usage and computational complexity by an order of magnitude. It can also converge the initial orbital error of hundreds of meters to the centimeter level within 10-15 iterations.

[0093] Overall architecture and lifecycle management of this invention (corresponding) Figure 1 , Figures 4-6 )

[0094] The solver adopts a four-level modular architecture: "foundation layer - algorithm layer - application layer - interface and management layer". All functions are encapsulated in a standard C++17 kernel, and C function interfaces prefixed with `orbit_solver_` are exported via `extern "C"`, forming a cross-language handle. During the initialization phase, the main process reads the initial orbital roots (...). Figure 4 ), station coordinates and inter-satellite / satellite-to-ground observation data ( Figure 5 After the spatiotemporal reference operator completes the ground-fixed to ground-inertial frame conversion, the data is injected into the base layer; simultaneously, a memory pool pre-allocation mechanism is utilized ( Figure 1 The base layer b) reserves the number of observations in advance to avoid heap fragmentation caused by dynamic expansion, achieving a data preparation stage with O(1) space complexity.

[0095] Distributed sparse matrices and field decomposition (corresponding) Figure 1 , Figure 2 , Figure 7 )

[0096] The base layer a uses Eigen::SparseMatrix<…> and extends support for the CSR-D dedicated format: row pointers only record satellite rows containing non-zero elements, and column indices are stored in segments consecutively according to the three perturbation sources: Earth's non-spherical shape, atmospheric drag, and solar radiation pressure. This ensures that items of the same type form continuous memory access in subsequent SPMVs, improving the L2 cache hit rate by ≥20%. Through "orbital plane-satellite" block row partitioning, each sub-matrix A_i corresponds one-to-one with the topology of domestic supercomputing nodes, with a fixed communication degree of 2 (…). Figure 2 This lays the foundation for overlapping communication.

[0097] RKF7(8) Adaptive Integral (corresponding to) Figure 11 )

[0098] The orbit calculation operator performs independent integration on a single satellite within each OpenMP thread, and the entire process is completed in the Earth inertial frame to avoid coordinate system confusion errors. The local truncation error Δ = ‖y7−y8‖∞ is estimated in real time. If Δ > ε_max, the step size is reduced by h_new = h·min(5, 0.8(ε_max / Δ)^(1 / 8)) to ensure that the single-satellite integration error is <1 cm. After the step size is accepted, the next step size is predicted so that the distribution of integration points adaptively matches the dynamic changes, reducing redundant calculations by about 25%.

[0099] Differential correction and Armijo line search (corresponding) Figure 3 , Figure 8 , Figure 12 )

[0100] After the normal equations are constructed, an Armijo inexact line search is performed along the direction d = −Δx. Figure 12 (Comparison chart). Traditional exact line search requires multiple global reductions to verify Wolfe conditions. This invention only retains one Iallreduce for the final step size α. The backtracking process is completely localized, and the number of communication synchronizations is reduced from 3 times per step to 1 time. In a 10,000-core scenario, the proportion of communication time is reduced from 42% to 12%, and CPU utilization is improved by ≥18%.

[0101] Mixed-precision PCG solution (corresponding) Figure 7 , Figure 13 )

[0102] The operator for solving linear systems automatically switches based on a sparsity threshold: PCG is enabled when the proportion of non-zero elements is ≥95%, the first 80% of iterations are performed using float, and the operator automatically upgrades to double when the relative change in residual δr < τ, with the final residual < 10⁻. 4 m, the computation time is reduced by ≥15%, the memory bandwidth usage is halved, and the numerical stability and computational efficiency are balanced.

[0103] Accuracy assessment and Allan variance sliding window (corresponding) Figure 10 )

[0104] The accuracy assessment operator uses the inverse matrix of the normal equation P=(AᵀWA)⁻¹ to propagate the position accuracy factor PPD=√tr(H·P·Hᵀ). An Allan variance sliding window is introduced (…). Figure 10 (Extended steps) σ_A(τ) is calculated in real time for the 30-epoch residual sequence, and τ_opt is automatically locked to realize dynamic PPD early warning, avoid sudden drop in accuracy, and improve the relative positioning accuracy of constellation by ≥20%.

[0105] Visualization and cross-language return (corresponding) Figure 1 Application layer c Figure 9 )

[0106] The visualization operators employ a single-line carriage return overlay technique in the command line, and the ProgressBar class dynamically refreshes the progress percentage in the terminal; the Statistics class calculates the mean and standard deviation online with a space complexity of O(1). orbit_solver_get_orbit() returns the orbit and PPD / VPD through a C structure, with a Python / Fortran call latency of <1ms, eliminating the need to recompile the upper-layer business code.

[0107] Fault-tolerant continuation calculation and exception handling (corresponding) Figure 1 Interface and management layer c)

[0108] When std::bad_alloc or matrix condition number > 1e12 is detected, the exception handling interface automatically triggers a checkpoint write-out; when resubmitting the job, orbit_solver_restore() reads the most recent checkpoint and resumes calculation within 30 seconds, improving the reliability of long-term tasks with tens of thousands of cores.

[0109] Furthermore, the distributed sparse matrix operation framework achieves communication hiding and storage optimization in the following ways:

[0110] The Jacobian matrix is ​​divided into blocks and rows according to the orbital plane-satellite number to form a submatrix A_i that corresponds one-to-one with the supercomputing node topology;

[0111] A double buffering mechanism is adopted. While performing sparse matrix-vector multiplication on the local submatrix A_i, the neighbor boundary vector Δx_{i±1} is exchanged through the MPI non-blocking interface Iallreduce, so as to achieve overlap of computation and communication.

[0112] Use MPI_Dist_graph_create_adjacent to build an adjacent communication graph of the orbital plane, fix the communication degree at 2, and reduce network congestion;

[0113] Based on the orbital dynamics Jacobi sparse structure, a dedicated CSR-D format is used for storage: row pointers only record satellite rows containing non-zero elements, and column indexes are stored in segments according to perturbation source type to reduce index memory and improve cache hit rate;

[0114] The analytically estimated block diagonal approximate inverse matrix M^{-1}_i is used as the preconditioner for domain decomposition to replace the dense Schur complement, thereby reducing the memory complexity of a single node to O(N).

[0115] Block partitioning and topology mapping

[0116] The Jacobian matrix of 10,000 satellites is divided into blocks in the order of "orbital plane-satellite number". Each block corresponds to the 6-dimensional state (position + velocity + dynamic parameters) of a satellite, resulting in a submatrix A_i. The partitioning rule corresponds one-to-one with the physical topology of the domestic supercomputing nodes: satellites in the same orbital plane are centrally assigned to the same MPI process group to ensure computation-data locality, avoid random access across nodes, and reduce the initial communication degree from O(N) to O(number of orbital planes).

[0117] Double-buffered computation-communication overlap

[0118] In each sparse matrix-vector multiplication (SPMV) stage, a "double buffering" mechanism is employed:

[0119] Buffer-A is used for the current computation core to multiply and accumulate local A_i·x_i;

[0120] Buffer-B is used to receive the neighbor boundary vector Δx_{i±1}.

[0121] The neighbor data exchange is initiated first through the MPI_Iallreduce non-blocking interface, and the Buffer-A calculation is executed immediately upon return. When the local calculation is completed, the Buffer-B has just finished receiving the data, and the latest boundary value can be used directly to update the residual. This achieves the overlap of calculation and communication time slices, and the global reduction delay is completely hidden. The measured communication time ratio is reduced by ≥40%.

[0122] Degree=2 Communication Graph Construction

[0123] By using MPI_Dist_graph_create_adjacent to add only adjacent orbital plane processes to the adjacency list, the in / out degree of each MPI process is fixed at 2, reducing the number of network links from O(P²) in the complete graph to O(P), significantly reducing the probability of switch congestion; at the same time, redundant handshakes are eliminated, and the Allreduce latency is weakly correlated with the number of processes, maintaining near-linear scaling even in scenarios with tens of thousands of cores.

[0124] CSR-D Dedicated Sparse Format

[0125] To address the "block sparse, intra-block dense" characteristics of the Jacobian in orbital dynamics, the CSR-D (CSR for Dynamics) format is designed:

[0126] The row pointer array only records satellite rows containing non-zero elements, skipping unobserved satellites, reducing index memory usage by approximately 35%.

[0127] The column index is stored in segments according to the perturbation source type (Earth is not spherical, atmospheric drag, solar radiation pressure), so that non-zero blocks of the same type are adjacent in memory. During SPMV, continuous loading is formed, the L2 cache hit rate is improved by ≥20%, and the single-node floating-point performance is improved by 18%.

[0128] Analyze block diagonal preconditioners

[0129] The analytically estimated block diagonal approximate inverse matrix M^{-1}_i=(A_ii)^{-1} is used as the preconditioner for domain decomposition, replacing the traditional dense Schur complement:

[0130] Storage capacity is linearly related to the number of satellites, and memory complexity has officially decreased from O(N²) to O(N).

[0131] When applying preconditioners, only local block inverse multiplication is required, eliminating the need for cross-node communication. The number of PCG iterations and the number of conditions remain at the same level as the original dense preconditioners, while the communication per iteration is reduced to zero, resulting in an overall solution time reduction of ≥30%.

[0132] By combining the above five steps, this invention can be realized on a domestically produced supercomputer with tens of thousands of cores:

[0133] Constant communication degree, hidden global reduction latency, and ≥78% efficiency for scaling to 10,000 cores;

[0134] The sparse index and cache-friendly design reduce the peak memory usage of the 8000-star task from the traditional 2.3 TB to 87 GB;

[0135] Preconditioner storage and computational complexity are linearized, supporting smooth scaling to constellations with tens of thousands of stars;

[0136] The overall orbit determination time is reduced by about 30%, providing a scalable, low-latency, and highly reliable linear algebraic foundation for the real-time operation and control of giant constellations.

[0137] Furthermore, when employing the conjugate gradient method, the linear system solver uses incomplete intra-node LU decomposition and inter-node domain decomposition preconditioners, and supports mixed-precision iteration, dynamically improving accuracy in the later stages of iteration to ensure numerical stability.

[0138] Furthermore, the visualization operator dynamically refreshes the progress percentage within a single line on the terminal using the ProgressBar class provided by progress_bar.hpp, and combines it with the Statistics class to provide real-time statistical output of satellite drag coefficient, orbital altitude, and residual RMS.

[0139] Furthermore, the management interface dynamically adjusts the load of the MPI process based on the satellite observation density. When the observation data of a certain orbital plane suddenly increases, the main process re-divides the satellite list and issues a new domain decomposition scheme to achieve dynamic load balancing.

[0140] Furthermore, the exception handling interface automatically triggers a checkpoint write when it detects a std::bad_alloc or matrix condition number exceeding the limit exception, and resumes iteration from the most recent checkpoint through the orbit_solver_restore() function when the task is resubmitted.

[0141] Mixed-precision PCG and preconditioner (corresponding) Figure 7 , Figure 13 )

[0142] Within a node: Incomplete LU decomposition (ILU(0)) is performed on a single A_i block, preserving the sparse pattern, with a fill factor < 1.2 and negligible memory increment; Between nodes: A domain decomposition preconditioner is used to approximate the Schur complement with the analytically estimated block diagonal inverse M⁻¹_i=(A_ii)⁻¹, avoiding cross-node communication. For the first 80% of the iterations, single-precision float is used, automatically upgrading to double when the relative change in residual δr < τ, with the final residual < 10⁻ 4 m, the computation time is reduced by ≥15% while ensuring numerical stability.

[0143] Visual operators (corresponding) Figure 9 , Figure 10 )

[0144] The ProgressBar class uses a single-line carriage return ('\r') to refresh the percentage and ETA in real time on the terminal; the Statistics class only maintains the online recursive formulas for the mean and variance, with a space complexity of O(1), and can output the mean and standard deviation of satellite drag coefficients, orbital altitude, and residual RMS in real time. It does not require storing all the data and is suitable for monitoring massive amounts of data at tens of thousands of stars.

[0145] Dynamic load balancing (corresponding to) Figure 2 , Figure 7 )

[0146] The main process periodically collects observation counts for each orbital plane. When the observation density of a certain plane is detected to be greater than the mean + 2σ, load balancing is triggered: the satellite list is re-partitioned, some satellites on the overloaded plane are moved to the lightly loaded nodes, and a new domain decomposition scheme is issued through MPI_Bcast. The computation and communication overlap, the wall-clock time of the balancing process is less than 5 seconds, and the CPU utilization rate is increased from 68% to 85%.

[0147] Exception handling and checkpoint recalculation (corresponding) Figure 1 Interface and management layer c)

[0148] The monitoring thread checks memory allocation and matrix condition count every 10 minutes. If std::bad_alloc is thrown or the condition count is greater than 1e12, orbit_solver_checkpoint_write() is automatically called to write the current state vector, covariance, and iteration count to the parallel file system. When the task is re-released, orbit_solver_restore() reads the most recent checkpoint and resumes the calculation within 30 seconds, avoiding starting a 10,000-core job from scratch and significantly improving the reliability of long-running tasks.

[0149] In another embodiment, a method for determining the orbit of a giant constellation using a decoder includes:

[0150] S1. Read the initial orbital elements, station coordinates, and inter-satellite / satellite-to-ground observation data;

[0151] This step specifically involves: file discovery and consistency verification.

[0152] The main process scans the specified directory at startup and identifies three types of input files based on predefined naming rules:

[0153] Initial orbital elements file (*.oe): Stores the Keplerian six elements and epoch time for each satellite;

[0154] Station coordinate files (*.stx): Store the Cartesian coordinates of all ground stations in the Earth-fixed system and the antenna height;

[0155] Observation data files (*.obs): Store inter-satellite / satellite-to-ground pseudoranges, pseudorange rate, carrier phase, and signal strength in a “satellite-station-epoch” hierarchy.

[0156] Integrity is verified by comparing the file magic number with the SHA-256 digest. If the data is missing or the hash is inconsistent, an ORBIT_SOLVER_DATA_CRC error is immediately thrown to avoid implicit deviations in subsequent calculations.

[0157] Parallel metadata parsing

[0158] MPI-IO is used to open all files collectively. The main process reads the header to obtain the total number of records N_sat, N_station, and N_obs, and then broadcasts this information to all processes by calling MPI_Bcast. Each process calculates the range of records [start, end) to be parsed locally according to the field decomposition scheme, avoiding redundant I / O. Zero-copy mapping is used in the parsing phase: file blocks are mapped to the virtual address space via mmap, and C++17 std::string_view slices each field to reduce one memory copy.

[0159] Unified conversion of spatiotemporal reference

[0160] The orbital elements obtained from the analysis are immediately fed into the spacetime reference operator:

[0161] Convert the epoch time from UTC to TAI, and then to relative seconds;

[0162] The station coordinates were transformed from ITRF2014 to the current epoch Earth-Fixed system using the IERS 2010 precession-nutation model, and then transferred to the Earth-Inertial system via the polar motion matrix.

[0163] For inter-satellite two-way ranging, construct Δρ_ij = ρ_ij − ρ_ji, eliminate common clock bias, generate "virtual observation" records, and provide symmetric row blocks for subsequent sparse matrix construction.

[0164] Memory pool pre-allocation and object reuse

[0165] The base layer memory pool management operator calls vector.reserve() to pre-allocate the orbit sequence and observation vector capacity once based on N_sat and N_obs obtained in step 2; at the same time, it reuses the stack buffer in the interpolation calculation to avoid frequent new / delete, so that the number of heap allocations in the 10,000-star data loading stage is reduced to constant level, and the loading time is shortened by ≥30%.

[0166] Data validity and gross error removal

[0167] The pseudorange and pseudorange rate were statistically tested using the 3σ criterion, and observation records with >3σ were marked as BAD_OBS. For satellites with more than 10 consecutive missing epochs, an ORBIT_SOLVER_DATA_GAP warning was triggered, and arc segments were automatically segmented to ensure a robust start for subsequent orbit determination iterations.

[0168] Cross-language handle injection

[0169] Finally, all the spatiotemporally unified and gross error-removed data is encapsulated into a C structure OrbitSolverData, which is injected into the solver kernel via orbit_solver_load_data(handle,&data). External Python / Fortran only needs to pass the file path string to complete the data loading, achieving zero-copy cross-language access.

[0170] S2. The basic layer operator converts the data into a ground inertial system state vector and constructs a sparse Jacobian matrix;

[0171] This step specifically involves a "zero-copy" rotational transformation from the solid to the inertial plane.

[0172] The basic spatiotemporal reference operator first transforms the station coordinates and orbital elements resolved in stage S1 from the ITRF2014 Earth-Fixed System to the current epoch Earth-Inertial System. Unlike the traditional double-buffered scheme of "converting to ECEF → ECI and then copying," this invention directly constructs the rotation matrix view in-situ on the read-only buffer of the mmap mapping: using C++17 consteval to generate the IERS 2010 precession-nutation-polar motion product matrix template at compile time, and performing rotation on the 128-bit aligned coordinate block at runtime in a single-instruction-multiple-data (SIMD) manner. The rotation result is immediately written to a pre-allocated temporary vector on the stack, avoiding an intermediate copy and reducing memory bandwidth usage by approximately 18%.

[0173] State vector parameterization extension

[0174] The classic six-element model is extended to a combined state vector x = [a, e, i, Ω, ω, M, Cd, Cr, dΩ / dt] consisting of six elements and dynamic parameters, where Cd and Cr are drag and radiation pressure coefficients, and dΩ / dt is used for long-term precession compensation. This extension is completed immediately after rotation, and the parameter block and orbital elements are stored contiguously in memory, providing a "one-time traversal, full partial derivative" basis for subsequent Jacobian block differentiation, reducing the number of repeated orbital propagations by 25%.

[0175] Sparse pattern "construction-filling" two-stage method

[0176] A three-level bitmap pre-scanning method based on "orbital plane-satellite-observation type" was creatively proposed:

[0177] Phase 1: Use a 16-bit bitmap to quickly mark which satellite-station pairs have observations at which epochs. The bitmap is stored compactly in orbital plane order, with a scanning complexity of O(N_obs) and is resident in the L1 cache.

[0178] Phase Two: The Jacobian matrix block sparse pattern is immediately calculated based on the bitmap, and the positions of non-zero blocks are determined in one step, avoiding branch prediction failures caused by the traditional "propagation-judgment" method. This two-phase method reduces the construction time of the 10,000-star sparse pattern from minutes to seconds.

[0179] Analysis - Numerical Mixed Partial Derivative Calculation

[0180] For the non-spherical J2, J3, and J4 terms, ∂a_J2 / ∂x is directly written using analytical partial derivative formulas, avoiding finite difference methods. For terms that are difficult to analyze, such as atmospheric drag and solar radiation pressure, a central difference method is used, but differentiation is only performed on the external parameters (Cd, Cr), while the orbital elements still reuse analytical results, balancing accuracy and speed. This hybrid strategy reduces the time required to construct Jacobian lines for a single satellite by 40% while maintaining meter-level difference accuracy.

[0181] Distributed CSR-D instant generation

[0182] After determining the sparse pattern, each MPI process generates a local CSR-D structure for the submatrix A_i: the row pointer array only records the rows containing observations; the column indices are stored contiguously in segments according to "Earth's non-spherical shape → atmospheric drag → solar radiation pressure," with non-zero blocks within the same segment tightly arranged to form a continuous memory access pattern. Compared with traditional CSR, this format saves 35% of index memory and improves the SPMV cache hit rate by 20%, laying a bandwidth foundation for subsequent PCG solving.

[0183] Neighbor boundary mapping and communication graph establishment

[0184] Utilizing the natural one-dimensional circular topology of the orbital plane, the adjacency relationship of each process is hard-coded as "left-right" neighbors by calling MPI_Dist_graph_create_adjacent, with a constant communication degree of 2. The resulting adjacency list is constructed only once during the program's lifetime, and all subsequent neighbor vector exchanges reuse this graph, eliminating the overhead of dynamic topology discovery.

[0185] Memory pool reuse and object lifecycle control

[0186] CSR-D triples (row_ptr, col_ind, values) and temporary partial derivative vectors are allocated from the base layer memory pool, using a "reserve first, then placement-new" strategy to ensure no extra heap allocation within the same orbital iteration; at the end of the iteration, reset() is called to zero the memory pool pointer but not to release it, reducing the number of system calls, and reducing the number of malloc / free calls during the construction of the 10,000-star matrix to constant level.

[0187] Through the above seven steps, the S2 stage completes the distributed construction of the ground inertial system state vector transformation and the 10,000-star sparse Jacobian matrix within seconds, providing a sparse and bandwidth-friendly linear system foundation for the subsequent S3 "integral-correction" closed loop, while maintaining cross-node index consistency and memory controllability.

[0188] S3. The algorithm layer operator performs RKF7(8) integration, Armijo differential correction and PCG solution until the convergence criterion is met;

[0189] This step specifically involves: an adaptive RKF7(8) integrator "single-core single-star" pipeline.

[0190] Each OpenMP thread exclusively uses one satellite, synchronously calculating k1–k8 using 7th and 8th order embedding formulas in the Earth inertial frame. The local truncation error Δ = ||y7−y8||∞ can be obtained immediately. If Δ > ε_max, the step size is reduced by h_new = h·min(5, 0.8(ε_max / Δ)^(1 / 8)); otherwise, the amplification factor for the next step is estimated. To avoid global barriers caused by step size differences between threads, a "step size autonomy + boundary interpolation" strategy is invented: the step size is determined autonomously within the thread, and interpolation is only needed at the observation time, eliminating the global synchronization bottleneck of the traditional integrator. The total integration time of 10,000 satellites increases approximately linearly.

[0191] Real-time feedback of observation residuals through "side integration-side construction"

[0192] An embedded residual construction module is used in the integration pipeline. When the integration time covers the observation epoch, the theoretical observation value is immediately calculated using the current state. The difference between the theoretical and measured values ​​is used to obtain the residual vector r, which is then written to the circular buffer in place. This buffer reuses stack memory to avoid frequent new / delete operations. The residual construction and integration are completed in the same loop, reducing one complete orbit extrapolation and saving approximately 18% of the overall CPU time.

[0193] Armijo's non-precise line search "single synchronization" mechanism

[0194] Traditional backtracking requires executing a global Allreduce to verify the Wolfe condition during each α trial. This invention retains only the Armijo sufficient descent condition Φ(x+αd)≤Φ(x)+c1α∇Φᵀd, and decomposes the calculation of the descent amount into:

[0195] Local components: Intra-thread accumulation of A_i^T r_i and (A_i d_i)^T (A_i d_i);

[0196] Global component: Only one Iallreduce is needed to reduce the total decrease and total norm after the final α is determined.

[0197] This reduces the number of communications per line search from 3 to 1, and the proportion of communication wall time in the 10,000-core scenario is reduced from 42% to 12%, while the backtracking process is completely localized, improving CPU utilization by ≥18%.

[0198] Analytical-numerical hybrid gradient update

[0199] For the non-spherical J2 / J3 / J4 terms of the Earth, analytical partial derivatives are used and written into the gradient block in one go; for difficult analytical terms such as atmospheric drag and solar radiation pressure, second-order central difference is used but the derivative is only taken on the external parameters (Cd,Cr), and the analytical results are reused for the orbital element part, which reduces the gradient construction time by 40% while maintaining meter-level difference accuracy.

[0200] Hybrid Precision PCG "Residual Trigger" Upgrade Strategy

[0201] The preconditioner uses intra-node ILU(0) + inter-node block diagonal approximation of the inverse M⁻¹_i, maintaining sparse storage and communication. For the first 80% of the iterations, it is executed using float, monitoring the relative change in residual δr = ||r_k|| / ||r_0||; when δr < τ (τ can be configured to 0.1), it is automatically upgraded to double, resulting in a final residual < 10⁻. 4 The computation time is reduced by ≥15%, memory bandwidth usage is halved, and the upgrade process does not require restarting the Krylov subspace, thus maintaining numerical stability.

[0202] Convergence Criterion "Dual Threshold" Adaptive

[0203] The outer iteration uses a residual RMS improvement threshold: convergence is determined when the RMS decrease ratio is less than 1% for two consecutive iterations and the mixed precision has been upgraded to double; the inner PCG uses a relative residual ||r_k||_M⁻¹<10⁻³. The dual-threshold mechanism avoids over-iteration, reduces the average number of iterations by about 2, and further shortens the overall solution time.

[0204] "Real-time output" at the end of the production line

[0205] Once the convergence criterion is met, the state vector, covariance, and PPD / VPD are immediately written into the pre-allocated structure, and the result file is written in parallel via MPI-IO. At the same time, the online Statistics class outputs the final RMS and drag coefficient mean, and the ProgressBar class displays 100% with automatic line breaks. The entire process is copy-free and barrier-free, ensuring a graceful exit of the 10,000-core job.

[0206] Through the above seven-step design, the S3 stage realizes an integrated pipeline of "integration-correction-solution" on the scale of 10,000 cores of domestic supercomputers: the number of communication times is significantly reduced, and the accuracy and stability are improved at the same time, providing core algorithm support for real-time centimeter-level orbit determination of giant constellations.

[0207] S4. Application layer operators output centimeter-level orbital data, PPD / VPD accuracy indicators, and visual statistical reports;

[0208] This step specifically involves: covariance-state synchronization propagation.

[0209] Once the PCG iteration of S3 achieves dual-threshold convergence, the parameter covariance matrix P=(AᵀWA)⁻¹ is immediately calculated on the local node. Using the pre-cached transformation matrix H=∂(r,v) / ∂x, the position-velocity covariance P_x=H·P·Hᵀ is obtained through a single SIMD concatenation multiplication method, and then the position precision factor PPD=√tr(P_rr) and velocity precision factor VPD=√tr(P_vv) are extracted. This propagation and state update are completed within the same MPI process, avoiding additional global reduction. The PPD / VPD latency is less than 1 ms, allowing it to be used by the upper layer.

[0210] In-situ centimeterization of orbital elements

[0211] The converged Cartesian state vector x=[r;v] of the Earth's inertial frame is converted in situ into Keplerian six roots. The conversion process reuses the angular momentum and energy vectors calculated in the integration stage, avoiding repeated square root and division operations, and the numerical error is controlled at the millimeter level. At the same time, the dynamic parameters Cd and Cr are extended and written into the same structure, so that the complete orbit and physical parameters are exposed to the outside world by a single pointer, and subsequent extrapolation does not require searching the database again.

[0212] Parallel file output and zero-copy disk writing

[0213] Employing the MPI-IO collective write mode, all processes simultaneously write local results by decomposing offsets by domain: orbital elements, PPD, VPD, and residual RMS fall into a single HDF5 dataset in a block-continuous manner; utilizing HDF5 chunks and compression filters, the amount of data written to disk in a single write for 10,000 stars is compressed to 60% of the original size, and the write time is reduced from minutes to seconds. Furthermore, the filename and dataset path are directly exposed through C-API, allowing Python / Fortran clients to open and read the dataset with zero copying.

[0214] Online Statistics O(1) Statistical Flow

[0215] All processes push the local residual sum of squares, drag coefficients, and orbital height to the Statistics class at the moment of convergence. This class only maintains three recursive variables: count, mean, and variance, with a space complexity of O(1) and no need to save historical data. The global mean and standard deviation can be obtained by a single MPI_Reduce(scalar, SUM), and the terminal outputs a streaming summary in real time, such as "Cd=2.157±0.023, RMS=0.031 m". The memory usage is constant, making it suitable for monitoring massive amounts of data at tens of thousands of stars.

[0216] Command-line ProgressBar single-row refresh

[0217] The ProgressBar class uses the ANSI escape sequence '\r' to loop and overwrite within a single line, decoupling the refresh frequency from the number of iterations. Combined with the global progress obtained from MPI_Reduce, it updates the percentage, ETA, and current PPD once per second, avoiding the screen-filling effect of traditional multi-line logs. Users can intuitively grasp the status of 10,000-core jobs in the SSH terminal without the need for a graphical interface or third-party monitoring software.

[0218] Dynamic accuracy warning and visualization curve

[0219] The sliding window Allan variance assessment is embedded at the end of S4: σ_A(τ) is calculated in real time for the residual sequence of the most recent 30 epochs, and the τ_opt that minimizes σ_A is automatically locked, thereby generating a dynamic PPD(t) curve; this curve is printed directly below the statistical summary in ASCII line form, achieving the effect of "terminal as graph". Engineers can visually identify the accuracy drift trend and trigger extrapolation or reorientation strategies in advance.

[0220] Cross-language struct return once

[0221] All results are ultimately encapsulated into a C structure OrbitSolverResult, which is returned to the external caller in one go via orbit_solver_get_result(handle, &result). The structure field order and memory alignment follow the ISO_C_BINDING specification. Fortran can directly use iso_c_binding for integration, while Python uses ctypes.Structure zero-copy mapping to eliminate the overhead of secondary serialization, reducing project integration time from "days" to "hours".

[0222] Through the above seven steps, the S4 stage completes centimeter-level orbit output, precision factor calculation, real-time statistics and visualization report generation within seconds. The entire process is without additional copying or global barriers, providing instant, accurate and easy-to-read orbit determination results for the online operation and control of giant constellations on domestic supercomputers.

[0223] S5. The interface layer returns track results and error codes via C functions, and supports external Python / Fortran calls.

[0224] This step specifically involves: handle lifecycle management.

[0225] The solver kernel is implemented in C++17, exposing only an opaque pointer to the outside world: `typedef struct orbit_solver_t* orbit_solver_handle;`. After instantiating the kernel object on the heap, `orbit_solver_create()` casts the original pointer to a handle and returns it. External callers only save the handle and cannot directly access the internal members, achieving "type erasure" and ABI isolation, avoiding cross-compiler incompatibility caused by C++ name mangling.

[0226] extern "C" exported function set

[0227] All external functions are declared with `extern "C"` to ensure stable symbol names. The core entry points include:

[0228] orbit_solver_load_data(handle, const char* path)

[0229] orbit_solver_run(handle)

[0230] orbit_solver_get_result(handle, OrbitSolverResult* out)

[0231] orbit_solver_get_last_error(handle, char* buf, size_t len)

[0232] Function parameters use only basic C types (pointers, size_t, double), eliminating STL, exceptions, and RTTI, ensuring direct linking for legacy languages ​​such as Fortran / Pascal.

[0233] Zero-copy result structure design

[0234] Define a C structure `OrbitSolverResult`: fields are arranged in descending order of memory alignment, with the first field being a `uint32_t magic` used for end-order verification. The track array, PPD, VPD, and RMS all use `double*` pointers to internal contiguous buffers; external callers obtain a "view" rather than a copy, avoiding secondary memory copying. The structure is appended with `size_tcount` and `size_t stride` to support Fortran array descriptors, implementing language-level zero-copy mapping.

[0235] Error Code - Context Chain Record

[0236] After a kernel exception is caught, it is converted into an OrbitSolverStatus enumeration code and written to a thread-local buffer containing the file line number, function name, and brief description. `get_last_error()` formats this information as `[ERR-0123] file:orbit.cpp:427, bad_alloc while allocating 8051MB` and returns, supporting location to the source file line and significantly shortening the debugging cycle for 10,000-core jobs.

[0237] Checkpoint and continuation handle serialization

[0238] `orbit_solver_checkpoint_write(handle, const char* path)` writes the current state vector, covariance, iteration count, and random number seed to a parallel HDF5 file. The path supports collective opening via `MPI_File_open`, ensuring contention-free simultaneous writing by thousands of processes. Upon resubmission, `orbit_solver_restore(handle, path)` reads the HDF5 file and reconstructs the C++ object, resuming computation within 30 seconds, achieving a "job-level snapshot".

[0239] Automatic generation of cross-language bindings

[0240] It provides the orbit_solver.h header file and the orbit_solver.modFortran module file, the latter of which is automatically generated by iso_c_binding, maintaining a 1:1 mapping between field names and C structures; on the Python side, it loads the .so file through ctypes.CDLL, and uses the pyorbit_solver.py wrapper to convert double* views to numpy.ndarray, enabling direct scientific computing without data copying, reducing the integration workload from days to hours.

[0241] Symbol visibility and version control: GCC__attribute__((visibility("default"))) is used to export only the orbit_solver_* prefix symbols, hiding the rest of the internal functions; at the same time, the ORBIT_SOLVER_VERSION macro and the runtime orbit_solver_version() string are provided to ensure that external callers can perform ABI version verification and avoid the risk of binary incompatibility during upgrades.

[0242] Through the above seven steps, the S5 interface layer achieves a complete solution of "one-time compilation, zero-copy multi-language calls, traceable errors, and resumable job calculations", enabling the giant constellation orbit determination kernel to be seamlessly embedded into legacy operation and control systems such as Beidou and Qianfan, significantly reducing engineering integration and maintenance costs.

[0243] Furthermore, the specific operational steps in step S3 are as follows:

[0244] a) Adaptive step size control: The local truncation error Δ = ||y is calculated using the RKF7(8) integrator. n +1(7)-y n +1(8‖∞,When Δ>ε_max, the step size is reduced by h_new=h·min(5, 0.8·(ε_max / Δ)^(1 / 8));

[0245] b) Armijo Inaccurate Line Search: Along the Correction Direction d k =-Δx backtracks, iteratively updating α←β·α (β=0.5) until Φ(x) is satisfied. k +αd k )≤Φ(x k )+c1α∇Φ k ᵀd k This reduces the number of global synchronized Allreduce operations by 70%.

[0246] c) Mixed-precision PCG solution: The first 80% of iterations use single-precision float, and the last 20% of iterations switch to double. Precision upgrade is automatically triggered by the relative change in residuals, and the final residual is <10⁻. 4 m and the calculation time is reduced by ≥15%.

[0247] a) Adaptive step size control

[0248] Each OpenMP thread independently executes the RKF7(8) embedding formula for a single satellite in the local inertial frame: within the same function call, the 7th order solution y7 is calculated first, and then the 8th order solution y8 is calculated. The local truncation error is instantaneously obtained as Δ=‖y7−y8‖∞. To avoid global barriers caused by differences in step size between threads, this invention adopts a "step size autonomy" strategy: when Δ>ε_max, the thread locally reduces the step size by h_new=h·min(5,0.8(ε_max / Δ)^(1 / 8)); when Δ<ε_min, the upper limit of the amplification factor for the same formula is 5 times. The amplified or reduced step size is limited to the hard window [h_min,h_max] to ensure numerical stability. When the integration time reaches the observation epoch, the thread calculates the state at that moment locally through three Hermite interpolations without additional orbit extrapolation, reducing the overall number of integrations by about 25%, while keeping the single-satellite position error <1 cm.

[0249] b) Armijo Inaccurate Line Search

[0250] The backtracking is performed along the correction direction d=−Δx. The core innovation lies in the "single synchronization" communication model: the traditional Wolfe condition requires global reduction of the gradient inner product in each α trial, while this invention only retains the Armijo sufficient descent condition Φ(x+αd)≤Φ(x)+c1α∇Φᵀd, and splits the descent amount into two parts: the locally accumulated A_i^T r_i and the global total descent. Specifically, the local component is calculated and accumulated in the register in each MPI process. Only when the backtracking terminates and the final α is determined, Iallreduce is executed once to reduce the total descent amount and the total norm. This reduces the number of Allreduce operations per linear search step from 3 to 1, and the communication wall time ratio in the 10,000-core scenario is reduced from 42% to 12%. The backtracking loop is completely localized, improving CPU utilization by ≥18%, and the non-precise characteristics do not compromise the final convergence accuracy.

[0251] c) Solving for mixed-precision PCG

[0252] The preconditioner employs a dual strategy: intra-node ILU(0) and inter-node analytic block diagonal inverse M⁻¹_i, maintaining sparse storage and communication. For the first 80% of the iterations, single-precision float is used, and the relative change in residual δr = ||r_k|| / ||r_0| is monitored in real-time. Once δr < τ (configurable to 0.1), the solver automatically switches to double, and the subsequent 20% of iterations are completed in double precision, resulting in a final residual < 10⁻¹. 4 The upgrade process does not require rebuilding the Krylov subspace; it only requires promoting float vectors to double in situ, instantly halving memory bandwidth and reducing overall solution time by ≥15%. This strategy significantly reduces memory usage and cache pressure while ensuring numerical stability, enabling near-linear scaling of PCG solutions on domestic supercomputing platforms.

[0253] In another embodiment, a computer-readable storage medium stores a computer program that, when executed by a processor, implements the steps of the above-described method.

[0254] In another embodiment, the present invention provides a computer-readable storage medium. This medium can be a parallel file system (such as Lustre, GPFS), general-purpose SSD, disk array, or cloud object storage of a domestically produced supercomputer, on which an executable program is persistently stored in binary form. The program source code is cross-compiled using a C++17 compiler and an MPI library to generate a dynamic shared library liborbit_solver.so with an extern "C" interface. This library contains all the instructions for implementing the orbit determination method for giant constellations: a complete process from data reading, geostationary-inertial system conversion, distributed sparse matrix construction, RKF7(8) adaptive integration, Armijo single synchronization line search, mixed-precision PCG solution to centimeter-level orbit output. When the processor loads and executes the program, it automatically completes the initialization of ten thousand cores in parallel processing, domain decomposition, communication graph establishment, and checkpoint continuation calculation. Without user intervention, the real-time orbit determination function for giant constellations can be implemented on the storage medium side, and centimeter-level precision orbit and PPD / VPD indicators can be output.

[0255] In another embodiment, an electronic device includes a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the program to implement the steps of the method described above.

[0256] In another embodiment, the present invention provides an electronic device. The device includes:

[0257] The storage device can be node-local DDR4 memory, NVMe SSD, or parallel file system mounted disk, used to store the aforementioned computer programs and multi-level input / output data.

[0258] The processor refers to the CPU computing core of a domestically produced supercomputer (such as Shenwei 26010 or Phytium FT-2000+), and the program is loaded into each core through the MPI launcher.

[0259] The network interface uses domestically produced Infiniband or Ethernet adapters for cross-node neighbor vector exchange and global reduction.

[0260] The checkpoint controller, called by the software layer, automatically writes the current state to memory when a memory anomaly or matrix condition number exceeds the limit, enabling continued calculation within 30 seconds.

[0261] When the processor executes a program in memory, the electronic device does so automatically:

[0262] Cross-language handle creation and data injection;

[0263] Construction and field decomposition of distributed sparse Jacobian matrices;

[0264] RKF7(8) adaptive integration, Armijo single-synchronous-line search, and mixed-precision PCG iteration;

[0265] Centimeter-level orbit, PPD / VPD accuracy factor, and real-time statistical report output;

[0266] Checkpoint writing and exception recovery.

[0267] When the entire electronic device is running at a scale of 10,000 cores, the number of communication synchronizations is reduced by 70%, the memory complexity is reduced from O(N²) to O(N), and the initial orbital error of hundreds of meters can be converged to the centimeter level in 10–15 iterations, thus forming a dedicated computing device with the ability to determine the real-time precise orbit of a giant constellation.

[0268] By adopting the above-disclosed technical solution of this invention, the following beneficial effects are obtained:

[0269] Scalability of 10,000 cores: Communication synchronization times are reduced by 70%, the measured scalability of domestic supercomputers is ≥78%, and the wall clock time for the same calculation is shortened by about 30%.

[0270] Memory controllable: Sparse format and region decomposition reduce the complexity of a single node from O(N²) to O(N), and the peak memory compression of the 8000-star task is over 90%, with no heap fragmentation during long-term operation.

[0271] Accuracy and stability: Adaptive step size + hybrid accuracy PCG final residual <10⁻ 4 m, converges to the centimeter level in 10–15 iterations; checkpoints continue calculation for 30 seconds, significantly improving the reliability of long-running operations with tens of thousands of cores.

[0272] Easy to use in engineering: C-API handles + zero-copy structures, zero recompilation and integration of Python / Fortran, reducing the integration cycle from "days" to "hours", and enabling "no graphical interface" monitoring of real-time progress and PPD curves on the terminal.

[0273] Resource utilization: Dynamic load balancing increases CPU utilization from 68% to 85%, reduces Infiniband link congestion, enables more orbit determination tasks to be completed with the same power consumption, and significantly reduces the overall cost of operation and control of giant constellations.

[0274] The above description is only a preferred embodiment of the present invention. It should be noted that for those skilled in the art, several improvements and modifications can be made without departing from the principle of the present invention, and these improvements and modifications should also be considered within the scope of protection of the present invention.

Claims

1. A constellation orbit determination solver based on a domestically produced supercomputer, characterized in that, The solver is implemented with a standard C++ kernel and exports a C-style function interface through extern "C". The interface manages the entire lifecycle of the solver using a handle mode. The solver is divided into four modules: a basic layer, an algorithm layer, an application layer, and an interface and management layer. The base layer includes: a) It adopts Eigen::SparseMatrix<...> and supports distributed sparse matrix storage in CSR, CSC, and COO formats; b) A memory pool management operator based on a memory pool pre-allocation + object reuse strategy, wherein the operator calls vector.reserve() to pre-allocate capacity before the track sequence is constructed and reuses the stack buffer during the interpolation process; c) A non-spherical gravitational perturbation model of the Earth, which includes at least J2, J3, and J4 band harmonic terms, and reserves an interface for extending spherical harmonic coefficients; d) Spatiotemporal reference operator to realize coordinate transformation between Earth-fixed system and Earth-inertial system and conversion between absolute calendar time and relative second; The algorithm layer includes: a) The orbit calculation operator uses the RKF7(8) adaptive step-size integrator. The integration process is completed entirely in the Earth inertial system. The step-size control estimates the local error by the difference between the 7th and 8th order solutions. b) For linear system solution operators, for the normal equation AT·A·Δx=–AT·r, the operation is automatically switched based on the non-zero element proportion threshold: when the sparsity is ≥95%, the conjugate gradient method is used; otherwise, LU decomposition is used. c) Differential correction operator, integrating inexact line search based on Armijo conditions, with initial step size α of 1.0 and decreasing backtracking until the sufficient descent condition is met; The application layer includes: a) The orbit determination operator transforms the giant constellation orbit determination problem into a closed-loop iteration of "prediction-comparison-correction". The convergence criterion for the iteration is that the improvement of the residual RMS is less than a set threshold for two consecutive iterations. b) Accuracy assessment operator: The parameter covariance matrix is ​​calculated based on the inverse matrix of the normal equation and propagated to obtain the position accuracy factor PPD and velocity accuracy factor VPD. c) Visual operators, using command line single-line carriage return overlay technology to output progress bars and online statistical summaries of Statistics in real time. The Statistics class outputs the mean and standard deviation with O(1) space complexity. The interface and management layer include: a) The programming interface is encapsulated, and all core functions are provided in the form of C functions prefixed with orbit_solver_. Error codes are returned using the OrbitSolverStatus enumeration. b) Management interface adopts a hybrid parallel architecture of MPI+OpenMP. MPI distributes the satellite constellation to different computing nodes according to the orbital plane. OpenMP parallelizes the propagation of single-satellite orbits within the node and realizes parallel file reading and writing through MPI-IO. c) Exception handling interface, providing get_last_error() to obtain the error code, location of occurrence and possible cause, supports checkpoint setting and recovery mechanism, and enables the continuation of long-term task interruption; The solver uses a domain decomposition method to distribute the sparse Jacobian matrix formed by the joint orbit determination of tens of thousands of satellites to domestic supercomputing nodes with tens of thousands of cores, thereby reducing the overall memory usage and computational complexity by an order of magnitude. It can also converge the initial hundred-meter-level orbital error to the centimeter level within 10-15 iterations.

2. The solver according to claim 1, characterized in that, The distributed sparse matrix operation framework achieves communication hiding and storage optimization in the following ways: The Jacobian matrix is ​​divided into blocks and rows according to the orbital plane-satellite number to form a submatrix A_i that corresponds one-to-one with the supercomputing node topology; A double buffering mechanism is adopted. While performing sparse matrix-vector multiplication on the local submatrix A_i, the neighbor boundary vector Δx_{i±1} is exchanged through the MPI non-blocking interface Iallreduce, so as to achieve overlap of computation and communication. Use MPI_Dist_graph_create_adjacent to build an adjacent communication graph of the orbital plane, fix the communication degree at 2, and reduce network congestion; Based on the orbital dynamics Jacobi sparse structure, a dedicated CSR-D format is used for storage: row pointers only record satellite rows containing non-zero elements, and column indexes are stored in segments according to perturbation source type to reduce index memory and improve cache hit rate; The analytically estimated block diagonal approximate inverse matrix M^{-1}_i is used as the preconditioner for domain decomposition to replace the dense Schur complement, thereby reducing the memory complexity of a single node to O(N).

3. The solver according to claim 1, characterized in that, When employing the conjugate gradient method, the linear system solver uses incomplete intra-node LU decomposition and inter-node domain decomposition preconditioners, and supports mixed-precision iteration, dynamically improving accuracy in the later stages of iteration to ensure numerical stability.

4. The solver according to claim 1, characterized in that, The visualization operator dynamically refreshes the progress percentage in a single line on the terminal using the ProgressBar class provided by progress_bar.hpp, and combines it with the Statistics class to output real-time statistics on satellite drag coefficient, orbital altitude, and residual RMS.

5. The solver according to claim 1, characterized in that, The management interface dynamically adjusts the load of the MPI process based on the satellite observation density. When the observation data of a certain orbital plane suddenly increases, the main process re-divides the satellite list and issues a new domain decomposition scheme to achieve dynamic load balancing.

6. The solver according to claim 1, characterized in that, When the exception handling interface detects a std::bad_alloc or matrix condition number exceeding the limit exception, it automatically triggers a checkpoint write and resumes iteration from the most recent checkpoint using the orbit_solver_restore() function when the task is resubmitted.

7. A method for determining the orbit of a giant constellation using the solver described in any one of claims 1-6, comprising: S1. Read the initial orbital elements, station coordinates, and inter-satellite / satellite-to-ground observation data; S2. The basic layer operator converts the data into a ground inertial system state vector and constructs a sparse Jacobian matrix; S3. The algorithm layer operator performs RKF7(8) integration, Armijo differential correction and PCG solution until the convergence criterion is met; S4. Application layer operators output centimeter-level orbital data, PPD / VPD accuracy indicators, and visual statistical reports; S5. The interface layer returns track results and error codes via C functions, and supports external Python / Fortran calls.

8. The method for determining the orbit of a giant constellation implemented by the solver according to claim 7, characterized in that, The specific steps in step S3 are as follows: a) Adaptive step size control: The local truncation error Δ = ||y is calculated using the RKF7(8) integrator. n +1(7)-y n +1(8‖∞,When Δ>ε_max, the step size is reduced by h_new=h·min(5, 0.8·(ε_max / Δ)^(1 / 8)); b) Armijo Inaccurate Line Search: Along the Correction Direction d k =-Δx backtracks, iteratively updating α←β·α (β=0.5) until Φ(x) is satisfied. k +αd k )≤Φ(x k )+c1α∇Φ k ᵀd k This reduces the number of global synchronized Allreduce operations by 70%. c) Mixed-precision PCG solution: The first 80% of iterations use single-precision float, and the last 20% of iterations switch to double. Precision upgrade is automatically triggered by the relative change in residuals, and the final residual is <10⁻. 4 m and the calculation time is reduced by ≥15%.

9. A computer-readable storage medium having a computer program stored thereon, the program being executed by a processor to implement the steps of the method of claim 8.

10. An electronic device comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor, when executing the program, implements the steps of the method of claim 8.

Images