%%%%%%%%%%%%%%%%%%%%%%% file template.tex %%%%%%%%%%%%%%%%%%%%%%%%%
%
% This is a general template file for the LaTeX package SVJour3
% for Springer journals.          Springer Heidelberg 2010/09/16
%
% Copy it to a new file with a new name and use it as the basis
% for your article. Delete % signs as needed.
%
% This template includes a few options for different layouts and
% content for various journals. Please consult a previous issue of
% your journal as needed.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% First comes an example EPS file -- just ignore it and
% proceed on the \documentclass line
% your LaTeX will extract the file if required
\begin{filecontents*}{example.eps}
%!PS-Adobe-3.0 EPSF-3.0
%%BoundingBox: 19 19 221 221
%%CreationDate: Mon Sep 29 1997
%%Creator: programmed by hand (JK)
%%EndComments
gsave
newpath
  20 20 moveto
  20 220 lineto
  220 220 lineto
  220 20 lineto
closepath
2 setlinewidth
gsave
  .4 setgray fill
grestore
stroke
grestore
\end{filecontents*}
%
\RequirePackage{fix-cm}
%
%\documentclass{svjour3}                     % onecolumn (standard format)
%\documentclass[smallcondensed]{svjour3}     % onecolumn (ditto)
%\documentclass[smallextended]{svjour3}       % onecolumn (second format)
\documentclass[twocolumn]{svjour3}          % twocolumn
%
\smartqed  % flush right qed marks, e.g. at end of proof
%
\usepackage{graphicx}
%\usepackage[linesnumbered,boxed,ruled,commentsnumbered]{algorithm2e}
\usepackage{color,array,amsthm}
\usepackage{soul, color, xcolor}
\usepackage[colorlinks,
linkcolor=blue, %% 修改此处为你想要的颜色
anchorcolor=blue, %% 修改此处为你想要的颜色
citecolor=blue %% 修改此处为你想要的颜色，例如修改blue为red
]{hyperref}
\documentclass{article}
\usepackage{graphicx}
\usepackage{subfigure}
\usepackage{float}
\usepackage{ragged2e}
\usepackage{amsmath, amssymb, mathrsfs}
\usepackage{marvosym}
\newcommand{\bigcorrespond}{\textsuperscript{\scalebox{1.2}{\Letter}}}
\usepackage{bbding}
\usepackage{pifont}
\usepackage{algorithm}
\usepackage{algpseudocode} 
\usepackage{breqn}
%
\usepackage{mathptmx}      % use Times fonts if available on your TeX system
\DeclareMathAlphabet{\mathcal}{OMS}{cmsy}{m}{n}
\DeclareSymbolFont{largesymbols}{OMX}{cmex}{m}{n}
%
% insert here the call for the packages your document requires
\usepackage{latexsym}
% etc.
%
% please place your own definitions here and don't use \def but
\newcommand{\familydefault}{\rmdefault}
%
% Insert the name of "your journal" with
\journalname{Nonlinear Dynamics}
%
\begin{document}

\title{Physics-informed Neural Network with Hard Constraints based on Dual-Step Training Strategy for Periodic Orbits in the CR3BP%\thanks{Grants or other notes
%about the article that should go on the front page should be
%placed here. General acknowledgments should be placed at the end of the article.}
}
%\subtitle{Do you have a subtitle?\\ If so, write it here}

%\titlerunning{Short form of title}        % if too long for running head

\author{Yong Chen\and
        Ming Cui \and Ying-Jing Qian * \and Ye-Feng Cheng \and Wen-Xue Chen %etc.
}

%\authorrunning{Short form of author list} % if too long for running head

\institute{Y. Chen \and M. Cui \and Y.J. Qian (\Letter) \and Y.F. Cheng \and W.X. Chen \at
              School of Mathematical Statistics and Mechanics, Beijing University of Technology, Beijing 100024, People’s Republic of China \\
              Corresponding author: Ying-Jing Qian (\email{candiceqyj@163.com})          %  \\
%             \emph{Present address:} of F. Author  %  if needed
}
%\date{Received: date / Accepted: date}
% The correct dates will be entered by the editor


\maketitle

\begin{abstract}
This study focuses on determining periodic orbits in the Circular Restricted Three-Body Problem (CR3BP) by introducing a novel Physics-informed Neural Network framework with hard constraints based on  Dual-Step Training (HDT-PINN). The proposed approach integrates several key components: initial hard constraints, a dual-step training strategy, and an augmented parameter optimization. The initial hard constraints enforce consistency between the computed orbit and the target orbit for the initial conditions. The  dual-step training strategy, including the pre-training step and final training step, effectively avoids convergence to local minima and is crucial for the successful identification of target periodic orbits. Furthermore, the augmented parameter optimization incorporates the initial velocity $\bar{\dot{y}}_0$ and the orbital period $t_f$ as trainable parameters, significantly enhancing the accuracy of the solutions. The effectiveness and robustness of the method are demonstrated across four distinct families of periodic orbits: Lyapunov orbits around the $L_1$ and $L_2$ points, distant retrograde orbits (DROs), and halo orbits. Furthermore, we conducted a convergence analysis of the initial velocity $\bar{\dot{y} }_0$ under various levels of disturbance and determined the convergence region for the initial velocity error $\delta_1$ associated with the periodic orbit considered in the example. The method proposed in this study exhibits a significantly broader convergence region in initial velocity error and precise target periodic orbit.
\keywords{Physics-informed neural network  \and Dual-step  training strategy\and CR3BP \and Periodic orbit \and Parameter optimization}
% \PACS{PACS code1 \and PACS code2 \and more}
%\subclass{MSC code1 \and MSC code2 \and more}
\end{abstract}

\section{Introduction}
%\label{intro}
%Your text comes here. Separate text sections with
%\section{Section title}
\label{sec: 1}
%Text with citations \cite{RefB} and \cite{RefJ}.

The Circular Restricted Three-Body Problem (CR3BP) within the Earth-Moon system remains a cornerstone of orbital dynamics, with profound theoretical significance and critical engineering implications for deep-space mission design and trajectory planning. In recent years, the Physics-informed Neural Network (PINN) has emerged as a powerful tool for solving scientific problems, demonstrating strong nonlinear fitting and equation-solving capabilities. Nevertheless, the application of PINN to the complex dynamics, particularly within frameworks like the CR3BP, remains underexplored, representing a significant area requiring rigorous exploration.


%\fontsize{10.0pt}{1.3em}\selectfont
Research on the CR3BP dates back to Newton. The foundational work by Poincaré \cite{HKC2001} in modern nonlinear theory demonstrated the absence of a general analytical solution for the three-body problem. The first systematic investigation of periodic orbits was conducted within the CR3BP framework, establishing them as fundamental to understanding the system's dynamics \cite{CELLETTI2025}. Based on this, analytical methods have been extensively studied in the CR3BP. A significant contribution was made by Howell \cite{HKC1984}, who introduced a Differential Correction (DC) technique that leverages a third-order analytical solution as an initial guess. Halo and Lissajous orbits were initially studied through the application of perturbation theory \cite{QYY2017}. More precisely, Farquhar \cite{FK1973} and Richardson \cite{RDL1980} derived third-order and fourth-order analytical approximations for Halo orbits near collinear libration points using the Lindstedt-Poincaré (LP) method. Furthermore, Peng et al. \cite{PENG2017} mainly studied the continuity of periodic orbits in the Sun-Mercury elliptical constrained three-body problem. Qian et al. \cite{QYZ2018,QYY2018} used the polynomial expansion method to analyze the periodic motion around libration points in the planar CR3BP and conducted parameter stability analysis for the Planar Bicircular Restricted Four-Body Problem (PBCR4BP). 

Since the 1980s, DC techniques have been widely applied in the CR3BP, which improves the accuracy of initial states and thus enhances the numerical accuracy of periodic orbits. This has attracted significant scholarly attention to numerical methods for addressing the CR3BP. Yuika \cite{Y2022} employed variational methods combined with symmetry constraints to derive periodic orbits in the planar CR3BP. Canales \cite{C2023} proposed an interpolation-based computational technique to identify periodic orbits, including both Lyapunov orbits and distant retrograde orbits (DROs). More recently, Tagliaferri et al.\cite{F2024} computed three dimensional periodic orbits near libration points using numerical approaches. Yamaguchi \cite{YGI2025} introduced an optimization strategy for Jacobi integral values, enabling the transfer of multiple asteroids into stable periodic capture orbits. Using the parameterization method within the PBR4BP, Peng et al. \cite{PLH2025} approximated the stable and unstable manifolds for unstable orbits near the Earth–Moon triangular libration points, accounting for solar gravitational perturbation.
While conventional methods have been foundational in CR3BP analysis, their potential for further breakthroughs may be constrained by certain methodological aspects. Specifically, their common approach of numerical integration from a single set of initial conditions links solution fidelity directly to the accuracy of those inputs. This linkage implies that computational orbits can be susceptible to divergence if the initial state contains errors.

With the development of artificial intelligence, machine learning has made significant progress in solving Partial Differential Equations (PDEs) in recent years \cite{YG2020}-\cite{WANG2026}. For instance, Rubinsztejn et al. \cite{AR2020} employed neural networks for optimal control in celestial dynamical systems. Moreover, a notable contribution in this field is the PINN framework, formally proposed by Raissi et al. \cite{RPK2019}, building upon earlier work by Lagaris et al. \cite{LLF1998}. The PINN methodology was developed primarily to address two classes of problems: data-driven solution of differential equations and data-driven discovery of governing equations. As a result, the PINN approach has garnered increasing research attention and has been successfully applied to fields such as celestial dynamics \cite{XDX2024}-\cite{ZL2024}.  Schiassi et al. \cite{EDD2022} applied PINN to orbital transfer in two-body problems. D'Ambrosio et al. \cite{DAE2021} achieved the orbit transfer from the halo orbit of libration point $L_1$ to that libration point $L_2$ using the PINN.
%\fontsize{10.0pt}{1.3em}\selectfont
However, the standard PINN shares the same drawbacks as traditional methods and still relies on theprecise initial values that can generate the target periodic orbit. Moreover, due to the strong nonlinearity inherent in the CR3BP, local optima pose a major challenge during the training process of standard PINN.

To address these dual challenges—the dependence on precise initial values and the prevalence of local minima—we propose a novel Physics-informed Neural Network with Hard Constraints based on Dual-Step Training Strategy (HDT-PINN) framework. This framework is designed to alleviate the reliance on precise initial conditions and enhance robustness against the CR3BP's strong nonlinearity. Our core improvements include a dual-step training strategy (pre-training and final training), the imposition of hard constraints on initial conditions, and an integrated parameter optimization mechanism. We validate the method by computing various periodic orbits and analyze the convergence region concerning initial velocity errors.


The paper is organized as follows. In Section \ref{subsec: 2}, the CR3BP is introduced, along with conventional PINN-based solution approaches. In Section \ref{subsec: 3}, the dual-step training strategy, hard constraints, and parameter optimization of HDT-PINN are detailed, accompanied by a comprehensive presentation of the algorithmic workflow. In Section \ref{subsec: 4}, the performance of the proposed method is evaluated through numerical experiments on various types of orbital trajectories. In Section \ref{subsec: 5}, the key contributions and findings of this study are summarized.


\section{CR3BP and PINN} \label{subsec: 2}
In this section, we first introduce the CR3BP, then the basic knowledge about standard PINN in solving periodic orbits is elaborated.



\subsection{Circular Restricted Three-Body Problem}
\setcounter{equation}{0}


The aim of this work is to investigate the properties of motion in the Circular Restricted Three-Body Problem (CR3BP). The two primaries move on circular orbits with the same Kepler frequency around their common center of gravity, which is assumed to be fixed at the origin of the coordinates. The third body moves in the same plane under the gravitational field of the two primaries. The nondimensional masses of the two primaries are $(1-\mu)$ and $\mu$, where $\mu = \frac{m_2}{(m_1 + m_2)} $ is the mass ratio. $m_1 $ and $m_2$ mean the mass of the two primaries, respectively. As depicted in Figure \ref{santiwenti}, the rotating coordinate frame where the origin is at the gravity center of two primaries is defined such that:

The $x$-axis aligns with the barycenter-connecting vector of primaries, directed toward the smaller primary;

The $z$-axis parallels the orbital angular velocity vector of the primaries;

The $y$-axis completes the right-handed Cartesian frame.

Then, the total time-independent potential in the rotating coordinate frame is
\begin{equation}\label{1_2}
	U = \dfrac{(x^2+y^2)}{2}+\frac{1-\mu}{r_1}+\dfrac{\mu}{r_2}+\frac{\mu(1-\mu)}{2}
\end{equation}	
where 
\[r_1 = \sqrt{(x+\mu)^2+y^2+z^2}\]
\[r_2 = \sqrt{(x+\mu-1)^2+y^2+z^2}\]
are the distances to the respective primaries.
%\justifying
%\fontsize{10.0pt}{1.3em}\selectfont

%\fontsize{10.0pt}{1.3em}\selectfont
The equations describing the motion of the third body in the rotating coordinate frame
\begin{equation}\label{1_1}
	\left\{
	\begin{array}{lcl}
		\ddot{x}-2\dot{y}=\frac{\partial U}{\partial x}\\
		\ddot{y}+2\dot{x}=\frac{\partial U}{\partial y}\\
		\ddot{z}=\frac{\partial U}{\partial z}\\
	\end{array}
	\right
\end{equation}				
%\justifying
%\fontsize{10.0pt}{1.3em}\selectfont

Furthermore, the position of the five Lagrangian points $L_1,\dots,L_5$ and the centers of the two primaries are shown in Figure \ref{santiwenti}. Three of them, $L_1, L_2$, and $ L_3$ are collinear points located on the $x$ axis. $L_1$ is a local minimum of the potential. The stationary points $L_2$ and $L_3$ are saddle points. The points $L_4$ and $L_5$, on the other hand, are local maxima of the gravitational potential.
\begin{figure}
	\centerline{\includegraphics[width=18.5pc]{santiwenti}}
	\caption{Distributions of five libration points in the three-body system}
	\label{santiwenti}
\end{figure}



\subsection{Standard PINN}


%\justifying
%\fontsize{10.0pt}{1.3em}\selectfont 
Since PINN introduction, the PINN method has consistently demonstrated exceptional capability in solving various types of differential equations across scientific and engineering fields. The PINN for solving systems of Ordinary Differential Equations (ODEs) is shown in Figure \ref{Physics-informed neural network framework}. Based on this, the PINN is used to solve the CR3BP, Eq.(\ref{1_1}) and (\ref{1_2}). Eq.(\ref{1_1}) is first transformed into the ODE problem $F_i(t,X)=0, i=1,2,3$, where $X=(x,y,z,\dot{x},\dot{y},\dot{z})$ represents the state variables in the system, which is expressed as
\begin{equation}\label{eq3}
    %\fontsize{8.0pt}{1.1em}\selectfont
    %\centering
    \left\{
		\begin{array}{lcl}
			F_1(t,x,y,z,\dot{x},\dot{y},\dot{z}) = \ddot{x}-2\dot{y}-\frac{\partial U}{\partial x}\\
			F_2(t,x,y,z,\dot{x},\dot{y},\dot{z}) = \ddot{y} + 2\dot{y}-\frac{\partial U}{\partial y} \\
			F_3(t,x,y,z,\dot{x},\dot{y},\dot{z}) = \ddot{z}-\frac{\partial U}{\partial z}\\
		\end{array}
        \right
	\end{equation}
and the initial condition is given by	
\begin{equation}
	\begin{aligned}
		x(t_0) = x_0,y(t_0) = y_0, z(t_0) = z_0,\\
        \dot{x}(t_0) = \dot{x}_0,\dot{y}(t_0) = \dot{y}_0, \dot{z}(t_0) = \dot{z}_0\\
	\end{aligned}  
\end{equation}
where $t_0$ denotes the initial time.
\begin{figure}
	\centerline{\includegraphics[width=18.5pc]{PINN}}
	\caption{Physics-informed neural network framework}
    \label{Physics-informed neural network framework}
\end{figure}

In the PINN (as seen in Figure \ref{Physics-informed neural network framework}), the solution of Eq.(\ref{1_1}) is approximated by a Single Input and Multiple Outputs (SIMO), which takes temporal $t$ variables as inputs, the approximated solution $g_k^\theta(k = 1,\dots, 6)$ as outputs, $g_k (k = 1,\dots, 6)$ means the states $(x,y,z,\dot{x},\dot{y},\dot{z})$, and  $g_{k0} (k = 1,\dots, 6)$ means the states $(x_0,y_0,z_0,\dot{x}_0,\dot{y}_0,\dot{z}_0)$. $\frac{\partial g_k^\theta}{\partial t}$ is implemented by using automatic differentiation, which can be easily realized in the deep learning framework like PyTorch. $\theta$ denotes the neural network parameter including the weights and biases of neural networks, whose value is computed by minimizing the loss function.
\begin{equation}
	\mathcal{L}_{total} = \omega_1 \mathcal{L}_{ODEs}+\omega_2\mathcal{L}_{p}+\omega_3\mathcal{L}_{0}
\end{equation}	
\begin{equation}
	\mathcal{L}_{ODEs}=\frac{1}{M}\sum_{j=1}^{M}\sum_{i=1}^{N_1}|F_i(t_j,g_k^\theta)|^2 
\end{equation}	
\begin{equation}
	\mathcal{L}_{p}=\sum_{k=1}^{6}|g_k^\theta(t_f)-g_k(t_f)|^2
\end{equation}	
\begin{equation}
	\mathcal{L}_{0}=\sum_{k=1}^{6}|g_k^\theta(t_0)-g_{k0}|^2
\end{equation}	
where  $\mathcal{L}_{ODEs}$ is the residual loss for a differential equations system, $t_f$ means the orbital period, $M$ is the number of points uniformly selected from the interval $[0, t_f]$, and $N_1$ is the number of differential equations, Enforcing $g_k(t_f)=g_{k0} $, $t_0=0$, then, $\mathcal{L}_{p}$ and $\mathcal{L}_{0}$ are the supervised learning loss at $t_f$ and $t_0$. Furthermore, $\mathcal {L}_ {total} $ is the total loss function constructed based on state variables, used to optimize the parameters $\theta$ of the neural network,  $\omega_1, \omega_2$, and $\omega_3$ mean the weight parameters. The optimization problem in the training of the PINN can be expressed as
\begin{equation}
	\left\{\theta,\hat{t}_f\right\} =\text{arg}\mathop{\text{ min}}\limits_{\theta,t_f}\mathcal{L}_{total}
\end{equation}	
where $\hat{t}_f$  denotes the optimized value.

%\justifying
%\fontsize{10.0pt}{1.3em}\selectfont
Despite the great achievement of the PINN method, it may still has the following issues in the strongly nonlinear CR3BP:

$\bullet$ The loss function is more complex, and the training process may fall into the trap of the local minima, which leads to the accuracy of the results has decreased.

$\bullet$ In the training process, PINN tends to enforce solutions' smoothness by sacrificing the initial conditions accuracy. This regularization effect may become particularly pronounced in the CR3BP scenarios with periodic boundary constraints \cite{AAS2021,KA2022}.

$\bullet$ Due to the strongly nonlinear nature of the three-body problem itself and its sensitivity to the initial conditions, the cumulative error of multi-objective training may be increased. 

Therefore, the dual-step training strategy is proposed in this study, the conclusion of the full pre-training step within the time domain is followed by the subsequent initiation of the final training step. 




\section{PINN with Hard Constraints based on Dual-Step Training Strategy}\label{subsec: 3}
To overcome the inherent limitations of standard PINN — specifically their propensity for local minima in highly nonlinear regimes and the susceptibility of pre-training to error accumulation in nonlinear systems such as the CR3BP \cite{GYW2023}. The HDT-PINN framework is proposed to advance standard PINN through three key enhancements: the imposition of hard constraints on initial conditions, the integration of an dual-step training strategy including pre-training and final training, and the adoption of parameter optimization, as illustrated in Figure \ref{fig:The relationship between PINN, PT-PINN, and HDT-PINN}.

%\fontsize{10.0pt}{1.3em}\selectfont

$\bullet$ Dual-step training strategy integrating the pre-training step with the final training step for enhancing stability: The dual-step training strategy can alleviate local minimum problems during the solving process and enhance the ability of neural networks to handle strongly nonlinear dynamics systems.

$\bullet$ Hard constraints drive the consistency between initial conditions and target orbits: The standard PINN can only approximate an orbit relative to given initial values. In contrast,  embedding hard constraints directly into the network's architecture guarantees that the target orbit based on initial values is strictly satisfied. This approach therefore ensures full consistency between the computed trajectory and the target periodic orbit defined by the initial conditions.

$\bullet$ Parameter optimization for enhanced robustness and accuracy: Building upon the incorporation of hard constraints into the network's architecture, optimized parameters are introduced to adjust initial states containing errors, thereby enabling precise refinement of the initial conditions. This approach significantly reduces the dependency of periodic orbit solutions on  accuracy of initial states, while simultaneously achieving accurate determination of the target orbit corresponding to the initial conditions.


\subsection{Dual-step training strategy}
The framework of the proposed method is depicted in Figure \ref{fig:ecpt-pinn}. Furthermore, this dual-step training strategy comprises two primary computational steps:
\subsubsection{Pre-training step}

%\fontsize{10.0pt}{1.3em}\selectfont

\begin{figure*}
	\centering
	\includegraphics[width=0.9\linewidth]{HDT-PINN.png}
	\caption{PINN with Hard Constraints based on Dual-Step Training strategy }
	\label{fig:ecpt-pinn}
\end{figure*}
\begin{figure*}
	\centering
	\includegraphics[width=0.8\linewidth]{方向流程图}
	\caption{The relationship between PINN and HDT-PINN }
	\label{fig:The relationship between PINN, PT-PINN, and HDT-PINN}
\end{figure*}
%\justifying
%\fontsize{10.0pt}{1.3em}\selectfont

According to the standard PINN framework, the state variables of the system are approximated by a Neural Network (NN), which can be a shallow or deep structure. Inspired by Schiassi et al. \cite{EDD2022}, we adopt SIMO neural network to approach state variables of the system.

In Figure \ref{fig:ecpt-pinn}, $t$ is the input and $y_k^\theta(t),(k = 1,\dots, 6)$ is the output of the neural network. Therefore, the solution satisfying the initial state through neural network approximation is expressed as 
\begin{equation}
	y_k^\theta(t) = \bar{ y}_{k0} + \frac{t}{t_f}g_k^\theta(t)
\end{equation}
%\justifying
%\fontsize{10.0pt}{1.3em}\selectfont
where $\bar{y}_ {k0},(k = 1,\dots,6)=x_0,y_0,z_0,\dot{x}_0,\bar{\dot{y}}_0,\dot{z}_0$ represents the initial state of the neural network's input. Since $\bar{\dot{y}}_0$ may be not accurate, here we add an optimization parameter $\varepsilon_1$, The neural network output $y_k^\theta(t)$ is corrected as
\begin{equation}
	y_k^\theta(t) = y_{k0} + \frac{t}{t_f}g_k^\theta(t)
\end{equation}
where $y_{k0},(k = 1,\dots,6) =x_0,y_0,z_0,\dot{x}_0,\bar{\dot{y}}_{0} + \varepsilon_1,\dot{z}_0$ is an more the accurate initial state achieving the target orbit. In the pre-training step, we set $\varepsilon_1 = 0$, $g_k^\theta(t)=\sigma(\omega t + b)$, $\sigma$ denotes the activation function, $\omega$ and $b$ denote the corresponding weights and biases. Besides, the derivative term $\frac{\partial y_k^\theta}{\partial t}$ appearing in the differential equation is obtained by using the automatic differential technique.

The pre-training step mainly divides the overall training time domain into several parts, that is, evenly dividing the entire time domain into $0=t_0<t_1 \dots<t_l=t_f $. Then, the initial pre-training time domain is set to $[t_0, t_1] $, and the training data is taken as $\tau_{l}$,

\[\tau_{1}=\left\{t_i | 0<t_i^1<t_1\right\}_{i=1}^{M_{1}}.\]
where $M_1$ means the number of training points on $[t_0, t_1] $.


The corresponding optimization problem is designed as
\begin{equation}
	\theta_1=\text{arg} \mathop{\text{min}}\limits_{\theta}\mathcal{L}_{total}(\tau_{1},\theta),
\end{equation}
where $\theta_1 $ is the parameter obtained by optimizing the loss function on $[t_0, t_1]$. 

Next, we will use the parameters and data obtained from each pre-training time domain for subsequent optimization. Concretely speaking, we take the pre-training time domain as $[t_0,t_{j}]$ and the training data as $ \tau_{j}$ into the next training step.
\[\tau_{j}=\left\{t_i | 0<t_i^j<t_j\right\}_{i=1}^{M_{j}}.\]



%\justifying
%\fontsize{10.0pt}{1.3em}\selectfont 
For ODEs, the setting of the general loss function only contains the initial conditions $\mathcal{L}_{0}$ and differential equations $\mathcal{L}_{ODEs}$. Since we transform the initial conditions into mandatory conditions, we only include $\mathcal{L}_{ODEs}$ and the supervised conditions obtained during the training process in our pre-training step.
\begin{equation}\label{3.8}
	\mathcal{L}_{total} = \omega_4 \mathcal{L}_{ODEs}(\theta,F_i,\tau_{j}) + \omega_5 \mathcal{L}_{sp}(\theta_{j-1},\tau_{sp})
\end{equation}
where $\omega_4$ and  $\omega_5$ are the weights of these two loss items, respectively. Furthermore, $\mathcal{L}_{sp}$ and $\mathcal{L}_{ODEs}$ can be written as
\begin{equation}
	\mathcal{L}_{ODEs}=\frac{1}{M_l}\sum_{k=1}^{M_j}\sum_{i=1}^{N_1}|F_i(t_{k}^{j})|^2 
\end{equation}
\begin{equation}\label{3.10}
	\mathcal{L}_{sp}=\frac{1}{M_{j-1}}\sum_{i=1}^{6}\sum_{k=1}^{M_{j-1}}|(y_i^\theta(t_{k}^{j-1}))-y^{\theta_{j-1}}_i(t_{k}^{j-1})|^2
\end{equation}

The labeled data sets for training the $\mathcal{L}_{ODEs}$ and $\mathcal{L}_{sp}$ are denoted by
\begin{equation}
	\tau_{j}=\left\{t_{i}^{j}|0<t_{i}^{j}<t_j\right\}_{i=1}^{i=M_j}
\end{equation}
\begin{equation}
	\tau_{sp}=\left\{(t_{i}^{j-1},y_k^\theta(t_{i}^{j-1}))|0<t_{i}^{j-1}<t_{j-1}\right\}_{i=1}^{i=M_{j-1}}
\end{equation}
where $M_j$ is the number of selected data points in $t \in [0, t_j]$, $\tau_{j}$ is the training data of the current pre-training step, $\tau_{sp}$ is the training data that we have obtained in the previous pre-training steps, which is utilized to the current training. Until training to $[t_0, t_l]$ using the training data, we further formulate the subsequent optimization problem.
\begin{equation}
	\theta_l   = \text{arg} \mathop{min}\limits_{\theta}\mathcal{L}_{total}(\theta_l,\theta_{l-1},\tau_{l},\tau_{sp})
\end{equation}

Finally, when $t_l \geqslant t_f$, we complete the pre-training step.

\subsubsection{Final training step}
%\justifying
%\fontsize{10.0pt}{1.3em}\selectfont


During the pre-training step, an approximate solution is obtained over the interval $[0, t_f]$. However,  this pre-training result alone fails to converge to the desired target orbit, as demonstrated in Sec.\ref{4.2.1}. To address this, a second-stage refinement—denoted as the final training step—is applied, using a refined loss function to meet the target performance criteria. The proposed loss function takes the following form
\begin{equation}\label{22}
	\mathcal{L}_{total} = \omega_6 \mathcal{L}_{ODEs}(\theta,F_i,\tau_l) + \omega_7 \mathcal{L}_{p}(\theta,t_f)
\end{equation}
\begin{equation}
	\mathcal{L}_{p} = \sum_{i}^{6}|y_i^\theta(t_f)-y_i(t_0)|
\end{equation}
%\justifying
%\fontsize{10.0pt}{1.3em}\selectfont 
where $\omega_6$ and  $\omega_7$ are the weights of these two loss items, respectively. 


To enhance optimization precision, we introduce a refined loss function $\mathcal{L}_{p}$, which expresses the problem's periodic constraints $y_k(t_f) = y_k(t_0),k=1,\dots,6$. Furthermore, we augment the iteration count during training to achieve convergence thresholds. Concurrently, the initial period estimate $\hat{t}_f$ is progressively refined throughout the optimization process. Employing the refined loss function enables global refinement of the periodic solution. Simultaneously, the initial parameter $\varepsilon_1$ is optimized. Thus, yielding the final optimization formulation:
\begin{equation}
	\left\{\theta,\varepsilon_1,\hat{t}_f\right\} =\text{arg}\mathop{\text{ min}}\limits_{\theta,t_f,\varepsilon_1}\mathcal{L}_{total}
\end{equation}
%\justifying
%\fontsize{10.0pt}{1.3em}\selectfont

After completing the above training process, the optimized parameters are substituted into the neural network and the target orbit will be obtained.

%\justifying
%\fontsize{10.0pt}{1.3em}\selectfont

\subsection{Hard constraints and parameter optimization}

Throughout our entire training process, in addition to dual-step training strategy,  there are also two key components: hard constraints and parameter optimization. PINN with hard constraints was proposed by LuLu et al. \cite{LP2021,CBE2025}, which strictly applies initial boundary conditions to the network structure by modifying the network output.


%\justifying
%\fontsize{10.0pt}{1.3em}\selectfont

As illustrated in Figure \ref{fig:ecpt-pinn}, the output \( g_k^\theta(t) \) from the neural network is incorporated with the initial condition to derive \( y_k^\theta(t) \). This formulation inherently enforces exact satisfaction of the initial condition's  target orbit throughout the solution process.

To further promote global optimality of the results, an adjustable parameter \( \varepsilon_1 \) is introduced into the initial condition, defined as:
\[
\varepsilon_1 = 0.5 \cdot y_{k0} \cdot \sigma_1(\bar{\theta}_1).
\]

Moreover, in the final training step, the training time domain is extended from $[0, t_f] $ to $[0, t_f+\varepsilon_2] $, where $\varepsilon_2$ is represented as
\[
\varepsilon_2 =  0.3 \cdot T \cdot \sigma_1(\bar{\theta}_2).
\]
where \( \sigma_1 \) is generally implemented as the \(\tanh\) function. $T$ represents the period time of the traditional DC solution. It is worth noting that $\varepsilon_1$ and $\varepsilon_2$ represent additional trainable parameters beyond the neural network, whose joint optimization is essential for accurately converging to the target orbit. Moreover, the parameters \( \bar{\theta}_1 \) and \( \bar{\theta}_2 \) are external to the neural network and are not updated during the pre-training step. Instead, they are optimized exclusively in the final training step. 





\begin{algorithm}
\caption{Framework of HDT-PINN}
\label{alg:algorithm1}
\begin{algorithmic}[1]
\State \textbf{Input:}
\State \quad A neural network $y^{\theta}_k(t) = y_{k}(t)$ satisfying $y^{\theta}_k(t_0) = y_{k0}$
\State \quad Initial state optimized parameter $\varepsilon_1 = 0 ,\varepsilon_2 = 0$
\State \quad A set of collocation points $\{t_i\}_{i=1}^{l} \in [t_0,t_f]$
\State \quad Training set $\tau_{i}=\{t_{j}^i \mid 0<t_{j}^i<t_l\}_{j=1}^{M_i}$
\State \quad $\tau_{sp}=\{(t_{j}^{i-1},y_k^\theta(t_{j}^{i-1})) \mid 0<t_{j}^{i-1}<t_{i-1}\}_{j=1}^{M_{i-1}}$
\State \quad Iteration number Niter of L-BFGS algorithm \cite{RPJC1995}

\Statex
\State \textbf{Pre-training step}
\State Neural network parameters initialization
\State Solve the standard pre-training optimization problem Eq.~(\ref{3.8}) in the first interval $(0, t_1]$
\State Output the neural network parameters

\For{$i = 1$ \textbf{to} $l$}
    \State Generate the supervised learning data set $\tau_{sp}$ in $[0,t_{i-1}]$ and training data set $ \tau_{i}$ in $[0,t_i]$
    \For{$n = 1$ \textbf{to} $N_{iter1}$}
        \State Minimize the loss function $\mathcal{L}_{total}(\theta_i,\theta_{i-1},\tau_{i},\tau_{sp})$ by using ~~~~L-BFGS optimizer:
        \State $\{\theta_i\} = \arg \min\limits_{\theta}\mathcal{L}_{total}(\theta_i,\theta_{i-1},\tau_{i},\tau_{sp})$
    \EndFor
    \State Refresh the neural network parameter $\theta_{i}$ as the initial guess of the next step
\EndFor

\Statex
\State \textbf{Final training step}
\For{$i = 1$ \textbf{to} $N_{iter2}$}
    \State Initial state optimization parameter  $\varepsilon_1 \in (-0.5 {\dot{y}}_0,0.5 {\dot{y}}_0)$ and period optimization parameter $\varepsilon_2 \in [-0.3T,0.3T]$
    \State Generate the supervised learning data set $\tau = [0,t_f]$
    \State Use the neural network parameter $\theta_{k}$ from the last training step optimize Eq.~(\ref{22})
    \State $\{\theta,\varepsilon_1,\hat{t}_f\} = \arg \min\limits_{\theta,t_f,\varepsilon_1}\mathcal{L}_{total}$
\EndFor
\end{algorithmic}
\end{algorithm}

	
	\section{Numerical examples}\label{subsec: 4}
	
%\justifying
%\fontsize{10.0pt}{1.3em}\selectfont
    
This section presents the application of the proposed method in the CR3BP. First, the implementation of the algorithm within the CR3BP framework is described. Subsequently, the hyperparameter selections and key configuration choices adopted for solving the CR3BP are specified. To validate the effectiveness of the method, numerical simulations of various types of periodic orbits—including both planar and three-dimensional cases—are presented. Finally, a convergence analysis is performed through numerical simulations under different perturbations to the initial conditions of the orbits.


\subsection{HDT-PINN for Period Orbit of CR3BP}

 The initial value of periodic orbits in  the CR3BP is given by
\begin{equation}
	(x_0,y_0,z_0,\dot{x}_0,\dot{y}_0,\dot{z}_0)
\end{equation}
where $\bar{\dot{y}} _0 =\delta_1 + \dot{y} _0 $ is the initial value after being perturbed. The $\delta_1$ is the perturbation, and $\dot{{y}}_0 $ is a traditional DC solution .
We provide an optimization parameter $\varepsilon_1 \in [-0.5 \dot {y}_0, 0.5 \dot{y}_0]$ to obtain the accuracy initial state, that is, the given initial value is rewritten as 
\begin{equation}
	(x_0,y_0,z_0,\dot{x}_0,\bar{\dot{y}}_0+\varepsilon_1,\dot{z}_0)
\end{equation}

Based on the equations of motion Eqs.(\ref{1_2}) and (\ref{1_1}), the periodic condition is written as
\begin{equation}\label{29}
    \begin{array}{lcl}
		x(t_f)=x_0,y(t_f)=y_0,z(t_f)=z_0\\
		\dot{x}(t_f)=\dot{x}_0,\dot{y}(t_f)=\bar{\dot{y}}_0+\varepsilon_1,\dot{z}(t_f)=\dot{z}_0	
    \end{array}
\end{equation}

The orbital period \( t_f \) is included in the input, subjecting to perturbation \(\delta_2 \in [-0.1 T,0.1  T]\). During the final training step, the parameter \( \varepsilon_2 \in [-0.3T,0.3T] \) is introduced as an additive correction term, such that  
\[
t_f \rightarrow t_f + \varepsilon_2,
\]  
to refine the period toward an accurate value.

According to the framework of the HDT-PINN, state variables approximation are derived through constraint expressions:
	\begin{equation}\label{29}
    \fontsize{8.0pt}{1.1em}\selectfont
    \centering
		\left\{
		\begin{array}{lcl}
			x(t) =\dfrac{t}{t_f} \sigma \beta_1 + x_0\\
			y(t) =\dfrac{t}{t_f} \sigma \beta_2 + y_0 \\
			z(t) = \dfrac{t}{t_f}\sigma \beta_3 + z_0\\
			\dot{x}(t) = \dfrac{t}{t_f}\sigma \beta_4 + \dot{x}_0\\
			\dot{y}(t) = \dfrac{t}{t_f}\sigma \beta_5 + \dot{y}_0+\varepsilon_1\\
			\dot{z}(t) = \dfrac{t}{t_f}\sigma \beta_6 + \dot{z}_0\\
			\ddot{x}(t) = \dfrac{t}{t_f}\dot{\sigma} \beta_4 + \dfrac{1}{t_f}\sigma \beta_4\\
			\ddot{y}(t) = \dfrac{t}{t_f}\dot{\sigma} \beta_5 + \dfrac{1}{t_f}\sigma \beta_5\\
			\ddot{z}(t) = \dfrac{t}{t_f}\dot{\sigma} \beta_6 + \dfrac{1}{t_f}\sigma \beta_6\\
		\end{array}
		\right
	\end{equation}
    
The activation function $\sigma$ is defined as a trainable linear combination:
\[
 \sigma(x) = a_1 \cdot \text{sigmoid}(x) + (1-a_1) \cdot \tanh(x)
\]
where the coefficients $a_1$ are adaptive parameter optimized during the training process. The gradient ${\sigma}$ is obtained using  the automatic differentiation technique \cite{JMLR2015}, $\beta_j(j=1,\dots, 6) $ represent weight and bias parameters.
	The loss function is composed of the dynamics equations, which is designed as：
	\begin{equation}\label{}
		\left\{
		\begin{array}{lcl}
			\mathcal{L}_{\ddot{x}} = \ddot{x} - 2\dot{y} - \frac{\partial U}{\partial x}\\
			\mathcal{L}_{\ddot{y}} = \ddot{y}+2\dot{x} - \frac{\partial U}{\partial y}\\
			\mathcal{L}_{\ddot{z}} = \ddot{z} - \frac{\partial U}{\partial z}
		\end{array}
		\right
	\end{equation}
   
	  However, the pre-training step optimizes the solution thr- ough time evolution, which confines the result to a specific time interval. The inability to determine precise orbital cycle times prevents the enforcement of periodic boundary conditions during the pre-training step. Therefore, the periodic condition is only used in the final training of the entire periodic result.
	\begin{equation}
		\begin{array}{lcl}
			\mathcal{L}_p = |x(t_f)-x_0|+|\dot{x}(t_f)-\dot{x}_0| \\
			~~~~~ + ~~|y(t_f)-y_0| + |\dot{y}(t_f)-(\bar{\dot{y}}_0+\varepsilon_1)| \\
            ~~~~~ + ~~|z(t_f)-z_0| + |\dot{z}(t_f)-\dot{z}_0|
		\end{array}                
	\end{equation}
    
     According to the particularity of the periodic orbit initial values, that is, the initial values are usually $(y = 0,\dot{x}(0)=0,\dot{z}(0)=0)$ , we also add the semi-periodic condition, that is
	\begin{equation}\label{eq32}
		\mathcal{L}_{0.5p} =|y(0.5t_f)|+|\dot{x}(0.5t_f)| + |\dot{z}(0.5t_f)|    
	\end{equation}

    The primary objective of integrating a pre-training step in the HDT-PINN framework is to obtain a preliminary solution approximation without stringent accuracy requirements. Then, we divide the pre-training into five steps, which are training on $[0,0.2t_f],[0,0.4t_f],[0,$ $0.6t_f],[0,0.8t_f],[0,t_f]$, and the number of training points in each step is $M_1 = 300,M_2 = 500,M_3 = 700,M_4 = 900,M_5 = 1100$.


	

In addition to the parameters given in the solving process mentioned above, other parameters in the entire algorithm also include an iteration count of $N_{\text{iter1}} = 500$. Regarding the weight setting of the loss function, the weights $\omega_4$ and $ \omega_5 $ of Eq. (\ref{3.8}) are set as $\omega_4=1, \omega_5=0 $ in $[0,0.2t_f]$. In subsequent training, supervised data is generated, and the weights are set as $\omega_4=0.6, \omega_5=0.4 $.
	
	Following pre-training step, the model enters its final training step. Then, based on Eq.(\ref{eq32}), the specified loss function is rewirtten as
	\begin{equation}
		\mathcal{L}_{total} = \omega_6 \mathcal{L}_{ODEs} + \omega_7 \mathcal{L}_{p} +\omega_{8} \mathcal{L}_{0.5p}
	\end{equation}
    

	%\justifying
	%\fontsize{10.0pt}{1.3em}\selectfont

	The model is trained over the full orbital period using $N = 500$ collocation points. The weights of the loss function are set as $\omega_6=0.6, \omega_7=0.2, \omega_{8}=0.2$. Moreover, an iteration count of $N_{\text{iter2}} = 5000$ is found to be sufficient for achieving convergence to the desired error tolerance across all periodic orbits considered in this study.

    Through the above settings, we can obtain the final optimized parameters：
    \[ \left\{\beta_1,\beta_2,\beta_3,\beta_4,\beta_5,\beta_6,\varepsilon_1,\hat{t}_f\right\} =\text{arg}\mathop{\text{ min}}\limits_{\beta_j,t_f,\varepsilon_1}\mathcal{L}_{total}\]

    Furthermore, in the full training process, the parameter configuration of the neural networks is introduced as follows.  Each state variable is represented using a single hidden layer neural network with 30 neurons. The L-BFGS optimizer is used to optimize the parameters in order to make the loss function $\text{mean}{|L|} =\sqrt {\mathcal{L}_ {total}}$  minmize.


    As for the PINN method, its parameter settings are consistent with those utilized in the final training step of the HDT-PINN method. It is noteworthy that the PINN approach implemented here is functionally equivalent to the HDT-PINN framework without the dual-step training strategy, parameters optimization, and hard constraints.

    In essence, the HDT-PINN function serves as a parametric solver for orbital dynamics. It takes an initial condition as input and, through a single forward pass that inherently encodes the physics‑informed training, directly outputs the complete orbit.



	%\justifying
	%\fontsize{10.0pt}{1.3em}\selectfont
	
	
	%\justifying
	%\fontsize{10.0pt}{1.3em}\selectfont
    

\begin{figure*}
\centering
\subfigure[Loss function of $L_1$ Lyapunov orbits]{\label{fig:Loss_example_1}
\includegraphics[width=0.45\linewidth]{Loss_example_1}}
\hspace{0.01\linewidth}
\subfigure[Loss function of $L_2$ Lyapunov orbit]{\label{fig:Loss_example_2}
\includegraphics[width=0.45\linewidth]{Loss_example_2}}
\vfill
\subfigure[$L_1$ Lyapunov orbit x-y]{\label{fig:example 2x-y1}
\includegraphics[width=0.45\linewidth]{example 1x-y_1}}
\hspace{0.01\linewidth}
\subfigure[$L_2$ Lyapunov orbit x-y]{\label{fig:example 1x-y1}
\includegraphics[width=0.45\linewidth]{example 2x-y_1}}
\caption{$L_1$ and $L_2$ Lyapunov orbit Loss function results and orbit solution results}
\label{L1 and L2 orbit error reduction results and orbit solution results}
\end{figure*}

    
    \begin{figure*}
\centering

\subfigure[Loss function of Low-Amplitude DRO]{\label{fig:Loss_example_3}
\includegraphics[width=0.45\linewidth]{Loss_example_3}}
\hspace{0.01\linewidth}
\subfigure[Loss function of Medium-Amplitude DRO]{\label{fig:Loss_example_4}
\includegraphics[width=0.45\linewidth]{Loss_example_4}}
\vfill
\subfigure[Low-Amplitude DRO x-y]{\label{example 3x-y_1}
\includegraphics[width=0.45\linewidth]{example 3x-y_1.png}}
\hspace{0.01\linewidth}
\subfigure[Medium-Amplitude DRO x-y]{\label{fig:example-4x-y1}
\includegraphics[width=0.44\linewidth]{example 4x-y_1.png}}
\caption{Low-Amplitude DRO and Medium-Amplitude DRO orbit Loss funcation results and orbit solution results}
\label{fig:Calculation results of planar periodic orbit}
\end{figure*} 


	\subsection{Periodic orbit results} \label{subsec: 4.2}
    
  
	%\justifying
	%\fontsize{10.0pt}{1.3em}\selectfont
	This section primarily discusses the applications of PINN, and HDT-PINN in various types of orbits, including $L_1$ Lyapunov Orbit, $L_2$ Lyapunov Orbit,  low-amplitude Distant Retrograde Orbit (DROs), medium-amplitude DROs, and Halo orbits. For each orbit type, a representative case is selected to analyze the convergence behavior of the loss function across different PINN-based architectures, the resulting orbital solutions, and the relative advantages and limitations of each PINN-based during the solution process.

        Tables \ref{L1 Lyapunov Orbit initial state}–\ref{Halo initial state} summarizes the initial conditions for the different orbit types. Both the standard PINN and HDT-PINN frameworks incorporate initial value optimization. To assess this capability, perturbed initial values—shown in the second row of each table—are introduced. Due to the fact that traditional orbital design methods usually employ the differential correction (DC) method to modify the analytical solution for obtaining initial values and subsequently converge to the target orbit. Thence, applying perturbation $\delta_1 $ to the velocity component $\dot{y}$ of the numerical solution after DC, a perturbation $\delta_2$ is applied to the orbital period $T$, as the input of the neural network in this study. Accordingly, an optimization parameter $\varepsilon_1$ is introduced to optimize $\bar{\dot{y}}_0$, and an optimization factor $\varepsilon_2$ is introduced to optimize $t_f$. The third and forth rows of each table present the optimized values of $\hat{\dot{y}}_0$ and $\hat{t}_f$ obtained via HDT-PINN and PINN methods, respectively.

     

	\begin{table*}
    \fontsize{8pt}{1.0em}\selectfont
		\caption{$L_1$ Lyapunov orbit initial states }
        \centering
		\label{L1 Lyapunov Orbit initial state}
		\begin{tabular}{ccccc}		
			\hline
			$L_1$ Lyapunov Orbit & $(x_0,y_0,\dot{x}_0,\dot{y}_0)$& $t_f$  & [\delta_1,\delta_2] \\
			\hline
			&  &  &    \\	
			Traditional DC Solution     &(0.8,0,0,0.35748)   & 3.3200 &  --\\
			&  &  &    \\	
			Initial Conditions with Error         & (0.8,0,0,0.3)    & 3.0  &    $[-0.16\dot{y}_0,0.096T]$ \\
			&  &  &  \\	
			HDT-PINN Optimized &(0.8,0,0.3573)   & 3.32001 &  -- \\
            &  &  &  \\	
			PINN Optimized &(0.8,0,0.3569)   & 3.3214 &   -- \\
			\hline
		\end{tabular}
	\end{table*}
    \begin{table*}
    \fontsize{8pt}{1.0em}\selectfont
		\caption{$L_2$ Lyapunov orbit initial states }
        \centering
		\label{L2 Lyapunov Orbit initial state}
		\begin{tabular}{ccccc}		
			\hline
			$L_2$ Lyapunov Orbit & $(x_0,y_0,\dot{x}_0,\dot{y}_0)$& $t_f$  & [\delta_1,\delta_2] \\
			\hline
			&  &  &    \\	
			Traditional DC Solution   &(1.15002,0,0,0.03011)   & 3.375 &  -- \\
			&  &  &    \\	
			Initial Conditions with Error       & (1.15002,0,0,0.025)    & 3.0    &  $[-0.169\dot{y}_0,0.11T]$   \\
			&  &  &    \\	
			HDT-PINN optimized &(1.15002,0,0,0.03014)   & 3.373 &    -- \\
            &  &  &    \\	
			PINN optimized &(1.15002,0,0,0.03011)   & 3.369 &    -- \\
			\hline
		\end{tabular}
	\end{table*}
    \begin{table*}
    \fontsize{8pt}{0.9em}\selectfont
		\caption{Low-Amplitude DRO initial states }
        \centering
		\begin{tabular}{ccccc}		
			\hline
			Low-Amplitude DRO  & $(x_0,y_0,\dot{x}_0,\dot{y}_0)$& $t_f$  & [\delta_1,\delta_2] \\
			\hline
			&  &  &   \\	
			 Traditional DC Solution  &(1.0406,0,0,-0.5373)   & 0.628 & -- \\
			&  &  &    \\	
			Initial Conditions with Error          & (1.0406,0,0,-0.51)    & 0.6    & $[-0.05\dot{y}_0,0.044T]$  \\
			&  &  &   \\	
			HDT-PINN Optimized&(1.0406,0,0,-0.5370)   & 0.629 & --   \\
            &  &  &   \\	
			PINN Optimized&  (1.068,0,0.25,-0.1943)   & 0.5 & --   \\
			\hline
		\end{tabular}
	\end{table*}

		\begin{table*}
        \justifying
	    \fontsize{8pt}{1.0em}\selectfont
		\caption{ Medium-Amplitude DRO initial states}
        \centering
		\begin{tabular}{ccccc}		
			\hline
			Medium-Amplitude DRO  & $(x_0,y_0,\dot{x}_0,\dot{y}_0)$& $t_f$ &[\delta_1,\delta_2]  \\
			\hline
			&  &     \\	
			 Traditional DC Solution    &(1.5609, 0, 0, -1.0505)   & 6.0958 & --  \\
			&  &      \\	
			Initial Conditions with Error         & (1.5609, 0, 0, -0.95)     & 6.0   & $[-0.09\dot{y}_0,0.01T]$   \\
			&  &     \\	
			HDT-PINN Optimized & (1.5609, 0, 0, -1.05049)   & 6.0944   &--  \\
            &  &     \\	
			PINN Optimized  & (1.5609, 0, -0.06, -1.0237)   & 6.1 & --  \\
			\hline
		\end{tabular}
	\end{table*}

    \begin{table*}
        \fontsize{8pt}{1.0em}\selectfont
   
		\caption{Halo orbit initial states }
        \centering
        \label{Halo initial state}
		\begin{tabular}{ccccc}		
			\hline
			Halo Orbit & $(x_0, y_0, z_0, \dot{x}_0, \dot{y}_0, \dot{z}_0)$& $t_f$  & $[\delta_1,\delta_2]$ & \\
			\hline
			&  &   &  \\	
			Traditional DC Solution    &(0.84966, 0, -0.17506, 0, 0.26306, 0)& 2.5598 &--\\
			&  &    &  \\	
			Initial Conditions with Error         & (0.84966, 0, -0.17506, 0, 0.22, 0)   & 2.5 & $[-0.163\dot{y}_0,0.022T]$\\
			&  &    & \\	
			HDT-PINN Optimized & (0.84966, 0, -0.17506, 0, 0.26303, 0)  & 2.5602 &--\\
            &  &    & \\	
			PINN Optimized     & (0.84966, 0, -0.17506, 0, 0.30774, 0)   & 2.995 &--\\
			\hline
		\end{tabular}
	\end{table*}

       
    \begin{figure*}
\centering

\subfigure[Loss function of Halo orbit]{\label{Loss_example_5.png}
\includegraphics[width=0.45\linewidth]{Loss_example_5.png}}
\hspace{0.01\linewidth}
\subfigure[ Halo orbit solution results]{\label{Halo_1.png}
\includegraphics[width=0.5\linewidth]{Halo.png}}
\caption{Halo orbit loss function error results and orbit solution results}
\end{figure*}  




	\subsubsection{$L_1$ Lyapunov Orbit and $L_2$ Lyapunov Orbit}\label{4.2.1}

Figures \ref{fig:Loss_example_1} and \ref{fig:Loss_example_2} exemplify the convergence behavior of the loss function when different PINN-based methods are applied to solve for $L_1$ and $L_2$ Lyapunov orbits. The convergence trends observed in these figures demonstrate the effectiveness of PINN-based approaches in computing periodic orbital solutions. Furthermore, HDT-PINN can achieve excellent results for optimization of initial states, according to the convergence of loss function.
    
     Figures~\ref{fig:example 2x-y1} and \ref{fig:example 1x-y1} present the trajectories of the $L_1$ and $L_2$ Lyapunov orbits computed using two methods: PINN and HDT-PINN methods initialized with perturbed states. When enhanced with parameter optimization, both the PINN and HDT-PINN methods are able to reconstruct target orbits with high precision. However, the proposed approaches demonstrate superior convergence and global exportation, as seen in the  Figures \ref{fig:Loss_example_1} and \ref{fig:Loss_example_2}.
    


   

	%\justifying
	%\fontsize{10.0pt}{1.3em}\selectfont

    
	

   

	

	

	\subsubsection{Low-Amplitude DRO and Medium-Amplitude DRO  }

	%\justifying
	%\fontsize{10.0pt}{1.3em}\selectfont
	
    
    Figures \ref{fig:Loss_example_3} and \ref{fig:Loss_example_4} illustrate the convergence behavior of the loss function across different PINN methods applied to DRO period orbit problems. The conventional PINN method tends to converge to local minima during optimization, failing to reach the globally optimal solution. In contrast, HDT-PINN methods employ dual-step training strategy, effectively avoid local minima and achieve convergence to the global optimum.

    Figures~\ref{example 3x-y_1} and \ref{fig:example-4x-y1} present the orbital trajectories those obtained from two PINN methods with perturbed initial conditions. The results demonstrate that the conventional PINN approach is prone to convergence to local minima during optimization, thereby failing to accurately reconstruct the target orbit. In contrast, the HDT-PINN method exhibits high numerical accuracy and robust stability in recovering both of the orbits under consideration.
    

    In summary, for the problem of solving for planar orbits, the standard PINN method can suffer from convergence to suboptimal solutions during optimization, limiting its applicability. In contrast, the proposed HDT-PINN method incorporates an enhanced dual-step training strategy coupled with parameter optimization. This approach effectively addresses the sensitivity of periodic orbits to initial conditions within the CR3BP and yields higher-fidelity numerical solutions for target orbits.



\FloatBarrier
 \begin{figure*}
\centerline{\includegraphics[width=33pc]{Figure_1111.png}}
\caption {Convergence region for different orbit types}
\label{Region of convergence range}
\end{figure*}
\FloatBarrier


  \begin{figure*}
\centering
\subfigure[$L_1$ initial value error orbit convergence result]{\label{fig:example 1}
\includegraphics[width=0.42\linewidth]{Figure_1.png}}
\hspace{0.01\linewidth}
\subfigure[$L_2$ initial value error orbit convergence result]{\label{fig:example 2}
\includegraphics[width=0.42\linewidth]{Figure_2.png}}
%\caption{$L_1$ and $L_2$ Lyapunov orbits convergence result with various initial value errors}
%\label{L1 and L2 orbit results and orbit solution results}

%\centering
\subfigure[$L_1$ initial value error orbit convergence result]{\label{fig:DC_1}
\includegraphics[width=0.42\linewidth]{DC_1.png}}
\hspace{0.01\linewidth}
\subfigure[$L_2$ initial value error orbit convergence result]{\label{fig:DC_2}
\includegraphics[width=0.45\linewidth]{DC_2.png}}
\caption{$L_1$ and $L_2$ Lyapunov orbits convergence result with various initial value errors}
\label{L1 and L2 orbit results and orbit solution results}
\end{figure*}

\FloatBarrier


  \begin{figure*}
\centering
\subfigure[Low-amplitude DRO initial value error convergence result]{\label{fig:example 3}
\includegraphics[width=0.42\linewidth]{Figure_3.png}}
\hspace{0.01\linewidth}
\subfigure[Medium-amplitude DRO initial value error convergence result]{\label{fig:example 4}
\includegraphics[width=0.42\linewidth]{Figure_4.png}}
%\caption{Low-amplitude DRO and medium-amplitude DRO orbit convergence result with various %initial value errors}
%\label{L1 and L2 orbit results and orbit solution results}
%\end{figure*}

%  \begin{figure*}
%\centering
\subfigure[Low-amplitude DRO initial value error convergence result]{\label{fig:DC_3}
\includegraphics[width=0.42\linewidth]{DC_3.png}}
\hspace{0.01\linewidth}
\subfigure[Medium-amplitude DRO initial value error convergence result]{\label{fig:DC_4}
\includegraphics[width=0.42\linewidth]{DC_4.png}}
\caption{Low-amplitude DRO and medium-amplitude DRO orbit convergence result with various initial value errors}
\label{DRO results and orbit solution results}
\end{figure*}

\FloatBarrier

   \begin{figure*}
\centering
\subfigure[Halo orbit  with various initial value errors]{\label{Halo 2.png}
\includegraphics[width=0.45\linewidth]{12345.png}}
\hspace{0.01\linewidth}
\subfigure[Loss function error of Halo orbit with various initial value errors]{\label{Loss_example 5.png}
\includegraphics[width=0.45\linewidth]{untitled1.png}}
\label{Halo orbit solution results}
\caption{Halo orbit convergence result with various initial value errors}

\end{figure*}  






	%\justifying
	%\fontsize{10.0pt}{1.3em}\selectfont
	
	
	\subsubsection{Halo Orbits}
    
	%\justifying
	%\fontsize{10.0pt}{1.3em}\selectfont


    
	In this subsection, the evaluation example of the proposed HDT-PINN algorithm is given in Halo orbits. Then, the initial states are shown in Table \ref{Halo initial state}. Figures \ref{Loss_example_5.png} and \ref{Halo_1.png} present the loss function convergence and resulting of orbits computed using the PINN and HDT-PINN. From the error convergence perspective, the standard PINN method remains prone to premature convergence to local minima. In contrast, the  HDT-PINN systematically address this issue for three-dimensional periodic orbits, thereby achieving and maintaining stable convergence. 
   

    As shown in Figure~\ref{Halo_1.png}, due to convergence limitations, traditional PINN cannot reconstruct the target orbit. In contrast, the proposed HDT-PINN method provides an algorithmic strategy that obtains more accurate target trajectories through perturbed initial numerical solutions after DC, while improving numerical stability throughout the entire training process.

   



	%\justifying
	%\fontsize{10.0pt}{1.3em}\selectfont






\subsection{Convergence region analysis}

To explore the convergence issues arising from inaccuracies in the initial conditions of the proposed method, we perform a systematic analysis of the convergence domains for different orbit types via detailed numerical simulations.  The initial conditions setup is same in Sec. \ref{subsec: 4.2}.  In this simulation analysis, we focus principally on the initial velocity. A uniform deviation of  $\delta_2 = -0.1T$ is applied to the initial period. For the initial velocity, the permissible error range is set to $[-0.5 \dot{y}_0, 0.5 \dot{y}_0]$. This configuration aligns with the optimization process, in which the velocity correction parameter $\varepsilon_1$ is constrained to the same normalized search interval, i.e., $\varepsilon_1 \in[-0.5 \dot{y}_0, 0.5 \dot{y}_0]$.

The convergence behavior under varying levels of perturbation in the initial velocity $\dot{y_0}$ is summarized in  Figure~\ref{Region of convergence range}. Furthermore, to visually and detailed illustrate the impact of varying initial perturbations on convergence, we present representative numerical results graphically, comparing cases where the perturbed initial conditions lie inside and outside the established convergence region, as shown in Figure~\ref{Region of convergence range}. 

The $L_1$ Lyapunov orbit corresponds to initial velocity errors of $-0.5\dot{y}_0$, $0.21\dot{y}_0$, and $0.22\dot{y}_0$, whereas the $L_2$ Lyapunov orbit is associated with errors of $-0.5\dot{y}_0$ and $0.5\dot{y}_0$, as illustrated in Figures~\ref{fig:example 1} and~\ref{fig:example 2}, respectively. Numerical results indicate that when the initial velocity error $\delta_1$ falls within the range of $[-0.5\dot{y}_0, 0.21\dot{y}_0]$, the target $L_1$ Lyapunov orbit can be successfully recovered. Notably, for the error interval $[0.22\dot{y}_0, 0.5\dot{y}_0]$, the solution converges to a distinct periodic orbit that is independent of the initial conditions, rather than restoring the predefined target orbit. In contrast, in Figure~\ref{fig:DC_1}, when $\delta_1=-0.5\dot{y}_0$ is applied, the DC method cannot retrieve the target orbit. In the L1 Lyapunov orbit, the correction range in the velocity error of traditional DC methods is significantly weaker than that of the HDT-PINN method.  For the $L_2$ Lyapunov orbit in the Figures~\ref{fig:DC_2}, successful recovery is achieved across the entire tested error range, i.e., for all $\delta_1 \in [-0.5\dot{y}_0, 0.5\dot{y}_0]$. This result is consistent with the predefined optimization parameter range $\varepsilon_1 \in [-0.5\dot{y}_0, 0.5\dot{y}_0]$, confirming that solutions for the $L_2$ orbit converge within the specified error bounds.

  

  The low-amplitude and medium-amplitude distant retrograde orbits (DROs) are illustrated in Figures~\ref{fig:example 3} and~\ref{fig:example 4}, respectively. Numerical results demonstrate that the target low-amplitude DRO can be successfully recovered when the initial velocity error $\delta_1$ lies within the range of $[-0.34\dot{y}_0, 0.5\dot{y}_0]$. For the medium-amplitude DRO, the target orbit is retrieved when $\delta_1$ falls within $[-0.09\dot{y}_0, 0.3\dot{y}_0]$. In contrast, for DROs, the traditional direct collocation (DC) method exhibits inferior performance to the HDT-PINN method under the condition of $\delta_1=-0.34\dot{y}_0$ (illustrated in Figure~\ref{fig:DC_3}). Nevertheless, the traditional DC method achieves excellent convergence performance for medium amplitude DROs, as shown in Figure~\ref{fig:DC_4}.

The corresponding calculated Halo orbits are presented in Figure~\ref{Halo 2.png}. Numerical findings indicate that the convergence region for the Halo orbit is defined as $\delta_1 \in [-0.5\dot{y}_0, 0.4\dot{y}_0]$. By comparison, the traditional DC method fails to converge to the target Halo orbit at $\delta_1 = -0.5\dot{y}_0$, and $0.4\dot{y}_0$, as seen in Figure \ref{Loss_example 5.png}. In the halo orbit, the convergence range of HDT-PINN is significantly higher than that of traditional DC methods, indicating that HDT-PINN has more stable performance and better convergence for the velocity error.
%\fontsize{10.0pt}{1.3em}\selectfont

Based on previous analysis of orbital convergence regions, the $L_2$ Lyapunov orbit exhibits the greatest convergence region of initial solution error. That is, given its stable convergence properties, it reliably converges to the target orbit under the fixed optimization parameter region  $\varepsilon_1 \in [-0.5 \dot{y}_0, 0.5 \dot{y}_0]$. In contrast, the medium-amplitude DRO demonstrates the narrowest convergence region, with stability only maintained for initial perturbations within $[-0.09 \dot{y}_0, 0.3 \dot{y}_0]$. Regardless of the extent of their respective convergence regions, these results demonstrate a crucial margin of robustness in the proposed method with respect to initial value constraints. This tolerance for imprecise or perturbed initial conditions underscores its practical utility where exact initialization is often infeasible. Furthermore, relative to traditional DC schemes, the presented HDT-PINN method exhibits a larger convergence region in the context of initial value perturbations, with the only exception being the case of moderately stable DRO. This finding confirms that the proposed method features enhanced universality and applicability under diverse orbital environments.

   \begin{table*}
    \fontsize{8.0pt}{1.0em}\selectfont
		\caption{The initial states of the $L_1$ Lyapunov orbit family.}
        \label{L1 Lyapunov orbit families}
		\centering
		\begin{tabular}{ccccccc}
	        \hline
			$(x_0,y_0,\dot{x}_0,\bar{\dot{y}}_0)$	&  $\text{mean}(|L|)$ & $\bar{\dot{y}}_0 + \varepsilon_1 $ & $\dot{y}_0$ & $t_f+\varepsilon_2$ & $T$ \\
	        \hline
            
	&  &  &  & & \\	
	(0.8000,0,0,0.3000)	& $9.0\times10^{-4}$   &0.3569& 0.357483 &3.32001 & 3.32\\
	&  &  &  &  \\	
	(0.805,0,0,0.28)&$9.0\times10^{-4}$    &0.3192& 0.319641 &3.124&3.142 \\
	&  &  &  &  \\	
	(0.81,0,0,0.25)	&$2.0\times10^{-3}$    &0.2714& 0.272288 & 2.980&3.000 \\
	&  &  &  &  \\	
	(0.815,0,0,0.2)	&$1.0\times10^{-3}$    &0.2174& 0.218012 &2.859&2.862 \\
	&  &  & &   \\	
	(0.82,0,0,0.15)	&$1.0\times10^{-3}$    &0.1624& 0.162513&2.780&2.780\\
	&  &  & &   \\	
	(0.825,0,0,0.1)	&$1.2\times10^{-3}$    &0.1093& 0.109819&2.730&2.750 \\
	&  &  &  &  \\	
	(0.83,0,0,0.05)	&$1.3\times10^{-3}$    &0.0610& 0.061125&2.703&2.705 \\
	\hline
\end{tabular}
\end{table*}

\begin{table*}
\fontsize{8pt}{1.0em}\selectfont
\caption{The initial states of the $L_2$ Lyapunov orbit family.}
\label{L2 Lyapunov orbit families}
\centering
\begin{tabular}{ccccccc}
         \hline
			$(x_0,y_0,\dot{x}_0,\bar{\dot{y}}_0)$	&  $\text{mean}(|L|)$& $\bar{\dot{y}}_0 + \varepsilon_1 $ & $\dot{y}_0$& $t_f+\varepsilon_2$& $T$ \\
	        \hline
&  &  &  &   \\	
(1.15478668, 0, 0, 0.0055)	& $2.0\times10^{-4}$   &0.00479& 0.004830 & 3.363&3.368\\
&  &  &   &  \\	
(1.15002190, 0, 0, 0.025)	& $3.0\times10^{-4}$   &0.03014& 0.030113& 3.373&3.375 \\
&  &  &   &  \\	
(1.145, 0, 0, 0.045)	& $4.0\times10^{-4}$   &0.05597& 0.056000& 3.377&3.374 \\
&  &  &  &   \\	
(1.14165858, 0, 0, 0.07)	& $5.0\times10^{-4}$   &0.07283& 0.072868& 3.379&3.377 \\
&  &  &   &  \\	
(1.13653079,0,0,0.09)	&$5.0\times10^{-4}$    &0.09823& 0.098303& 3.384&3.385\\
&  &  &   &  \\	
(1.13084989,0,0,0.11)	&$6.0\times10^{-4}$    &0.12584& 0.125975& 3.394&3.395\\
&  &  &   &  \\	
(1.12450898,0,0,0.13)	&$5.0\times10^{-4}$    &0.15642& 0.156410 & 3.406&3.4055\\
&  &  &  &   \\	
(1.11345329,0,0,0.2)	&$5.0\times10^{-4}$    &0.20888& 0.208877&3.444&3.445\\
&  &  &    \\	
(1.10927089,0,0,0.2)	&$2.0\times10^{-3}$    &0.22833& 0.228685&3.430&3.431\\
\hline
\end{tabular}
\end{table*}

    \begin{table*}
\fontsize{8.0pt}{1.0em}\selectfont
\caption{The initial states of the DRO family.}
\label{DRO families}
\centering
\begin{tabular}{ccccccc}
\hline
$(x_0,y_0,\dot{x}_0,\bar{\dot{y}}_0)$	&  $\text{mean}(|L|)$& $\bar{\dot{y}}_0 + \varepsilon_1 $ & $\dot{y}_0$& $t_f+\varepsilon_2$& $T$  \\
	        \hline
&  &  &    \\	
(1.35, 0, 0, -0.78)	& $8.0\times10^{-4}$   &-0.71296& -0.71292&5.358 &5.359 \\
&  &  &    \\	
(1.40, 0, 0, -0.85)	& $9.0\times10^{-4}$   &-0.78938& -0.78945 & 5.651& 5.650\\
&  &  &    \\	
(1.45, 0, 0, -0.9)	& $1.0\times10^{-3}$   &-0.86861& -0.86868& 5.849 & 5.849 \\
&  &  &  & & &  \\	
(1.50, 0, 0, -0.92)	& $1.0\times10^{-3}$   &-0.94983& -0.94980 & 5.983&5.984\\
&  &  &    \\	
(1.5609,0,0,-0.95)	&$1.2\times10^{-3}$    &-1.05054& -1.05057& 6.094&6.0958\\
&  &  &    \\	
(1.60,0,0,-1.1)	&$1.5\times10^{-3}$    &-1.11588& -1.11628 &6.144&6.146\\
&  &  &    \\	
(1.65,0,0,-1.15)	&$2.5\times10^{-3}$    &-1.20106& -1.20149& 6.191&6.193 \\
&  &  &    \\	
(1.70,0,0,-1.2)	&$3.0\times10^{-3}$    &-1.28813& -1.28823&6.224&6.227\\
\hline
\end{tabular}
\end{table*}



\begin{table*}
\fontsize{8.0pt}{0.9em}\selectfont
\caption{The initial states of the Halo family.}
\label{Halo families}
\centering

\begin{tabular}{ccccccc}
\hline
$(x_0,y_0,z_0,\dot{x}_0,\bar{\dot{y}}_0,\dot{z}_0)$	&  $\text{mean}(|L|)$& $\bar{\dot{y}}_0 + \varepsilon_1 $ & $\dot{y}_0$& $t_f+\varepsilon_2$& $T$\\
\hline
&  &  &    \\	
(0.86129,0,-0.18487,0,0.22,0)& $8.0\times10^{-4}$   & 0.25237& 0.25239 &2.381 &2.383 \\
&  &  &    \\	
(0.84966,0,-0.17506,0,0.22,0)	& $9.0\times10^{-4}$   &0.26303& 0.26306& 2.560& 2.558\\
&  &  &    \\	
(0.84332,0,-0.16525,0,0.22,0)	& $9.0\times10^{-4}$   &0.26373& 0.26376 &2.659 &2.662\\
&  &  &    \\	
(0.83926,0,-0.15544,0,0.22,0)	& $7.0\times10^{-4}$   &0.26031 & 0.26036& 2.716&2.713 \\
&  &  &    \\	
(0.8363,0,-0.1456,0,0.20,0)	&$1.0\times10^{-3}$    & 0.25458& 0.2549& 2.75&2.77\\
&  &  &    \\	
(0.8320,0,-0.1260,0,0.20,0)	&$4.0\times10^{-4}$    &0.24017& 0.2404& 2.782&2.780 \\
&  &  &    \\
(0.8303,0,-0.0162,0,0.20,0)	&$1.1\times10^{-3}$    &0.23148& 0.2314 & 2.787&2.790\\
&  &  &    \\
(0.8289,0,-0.1064,0,0.20,0)	&$5.0\times10^{-4}$    &0.22162& 0.2220&2.785 &2.783 \\
&  &  &    \\	
(0.8276,0,-0.0966,0,0.18,0)	&$9.0\times10^{-4}$    &0.21151& 0.2114&2.784 &2.783  \\
&  &  &    \\
(0.8264,0,-0.0868,0,0.18,0)	&$8.0\times10^{-4}$    &0.20127& 0.2013& 2.780&2.775  \\
&  &  &    \\	
(0.8255,0,-0.0770,0,0.16,0)	&$1.4\times10^{-3}$    &0.18990& 0.1903&2.774 &2.771\\
&  &  &    \\	
(0.8247,0,-0.0672,0,0.15,0)	&$9.0\times10^{-4}$    &0.17915& 0.1791& 2.768&2.77 \\
\hline
\end{tabular}
\end{table*}  
    
	\begin{figure*}
\centering

\subfigure[$L_1$ Lyapunov orbital family]{\label{fig:example-2x-y1}
\includegraphics[width=0.31\linewidth]{example 1}}
\hspace{0.01\linewidth}
\subfigure[$L_2$ Lyapunov orbital family]{\label{example 1x-y_1}
\includegraphics[width=0.31\linewidth]{example 2}}
\hspace{0.01\linewidth}
\subfigure[DRO family]{\label{example 3}
\includegraphics[width=0.31\linewidth]{example 3}}
\caption{Planar orbital family}
\label{Plane orbit families}
\end{figure*}
 \begin{figure*}
\centerline{\includegraphics[width=18pc]{example 4_2.png}}
\caption {Halo orbital family}
\label{fig:Distant retrograde orbit families}
\end{figure*}
   
	\subsection{Related orbital families}



	%\justifying
	%\fontsize{10.0pt}{1.3em}\selectfont
	To validate the generality of our method across various orbital types, we employed the HDT-PINN approach to compute representative members of the $L_1$ and $L_2$ Lyapunov orbital families, as well as the DRO orbital family and Halo orbital family. 

    Tables \ref{L1 Lyapunov orbit families}-\ref{Halo families} summarize the computational results for the \(L_1\) and \(L_2\) Lyapunov orbital families, DRO family, and Halo orbital family obtained via the HDT-PINN method. The first column of each table lists the selected initial state for the respective family. These initial conditions are derived from traditional DC method, with the initial velocity \(\dot{y}\) subjected to a perturbation \(\delta_1\). For the given \(\delta_1\), all trajectories successfully achieved the optimization objective of converging to the target orbit and refining the initial values. The column labeled "\(\text{mean} | L |\)" reports the training loss function for each orbit. The third column provides the optimized initial value \(\bar{\dot{y}} + \varepsilon_1\), and \(\dot{y}\) denotes the initial value from the reference traditional DC method. Here, \(\hat{t}_f = t_f + \varepsilon_2\) and \(T\) represent the optimized orbital period and the accurate orbital period, respectively.
    




    As our methodology incorporates periodic boundary conditions, the initial period \( t_f \) must also be treated as an optimization variable. In the orbit computations above, a perturbation \( \delta_2 \in [-0.2 T, 0.2 T] \) was applied to the orbital period \( T \) of the traditional DC solution . By introducing the optimization parameter \( \varepsilon_2 \), the complete periodic target orbit is successfully solved across all cases.



   During the optimization process, these parameters $\varepsilon_1$ and $\varepsilon_2$ are typically assigned relatively broad ranges, such as $\varepsilon_1 \in [-0.5 \dot{y}_0, 0.5 \dot{y}_0]$ and $\varepsilon_2 \in [-0.3T, 0.3T]$. This configuration enhances the convergence speed of the optimization and aids in the stable identification of the optimal initial velocity $\hat{\dot y}_0 = \bar{\dot{y}}_0 + \varepsilon_1$ and period $\hat{t}_f = t_f + \varepsilon_2$ . Finally, the resulting optimization accuracy is as follows:
$\left| \dot{y}_0 - (\bar{\dot{y}}_0 + \varepsilon_1) \right| < 8 \times 10^{-4}$ and $\left|T - (t_f + \varepsilon_2) \right|$$ < 5 \times 10^{-3}$.



    The results for the \( L_1 \) and \( L_2 \) Lyapunov, DRO, and Halo orbital families are visualized in Figures \ref{Plane orbit families} and \ref{fig:Distant retrograde orbit families}. These results demonstrate that the HDT-PINN method possesses the capability to solve for a diverse set of orbital families. Furthermore, it exhibits notable generalization and precision performance for orbital determination within the algorithm framework.








%		\begin{tabular}{ccccc}
%	\hline
%	$(x_0,y_0,\dot{x}_0,\dot{y}_0)$	&  $\text{mean}(|L|)$& $\dot{y}}_0 + %%\epsilon_1$&$T + \epsilon_2$& $\dot{\hat{y}}}_0$  \\
%		\hline
%		&  &  &  &  \\	
%		(0.8,0,0,0.3)	& $9.0\times10^{-4}$   &0.3569&3.3200& 0.357483 \\
%		&  &  &  &  \\	
%		(0.805,0,0,0.28)&$9.0\times10^{-4}$    &0.3192&3.1249& 0.319641 \\
%		&  &  &  &  \\	
%		(0.81,0,0,0.25)	&$2.0\times10^{-3}$    &0.2714&2.9806& 0.272288 \\
%		&  &  &  &  \\	
%		(0.815,0,0,0.2)	&$1.0\times10^{-3}$    &0.2174&2.8594& 0.218012 \\
%		&  &  &  &  \\	
%		(0.82,0,0,0.15)	&$1.0\times10^{-3}$    &0.1624&2.7801& 0.162513\\
%		&  &  &  &  \\	
%		(0.825,0,0,0.1)	&$1.2\times10^{-3}$    &0.1093&2.7300& 0.109819 \\
%		&  &  &  &  \\	
%		(0.83,0,0,0.05)	&$1.3\times10^{-3}$    &0.0610&2.7030& 0.061125 \\
%		\hline
%\end{tabular}




\section{Conclusions} \label{subsec: 5}
%\justifying
%\fontsize{10.0pt}{1.3em}\selectfont
In this study, we employ a PINN with hard constraints based on the dual-step training strategy, integrating pre-training step and final training step, to compute periodic orbits in the CR3BP. Specifically speaking, the dual-step training strategy is adopted to mitigate the issue of local minima that often hinders conventional PINN-based solvers. Furthermore, the hard constraints on the initial states and the augmented parameter optimization are introduced to address the high sensitivity of the CR3BP to initial conditions, which enhances the convergence ability and the overall accuracy of the target orbit corresponding to the initial conditions. By imposing periodic conditions as fundamental constraints, the proposed method achieves accurate reconstruction of the target orbit even under perturbed initial states. 


To evaluate the performance of the proposed approach, multiple orbit types are examined. The results demonstrate that the HDT-PINN method not only remains robust under initial conditions with perturbation in orbit determination, but also improves the estimation of the initial state. Therefore, this proposed method enabling successful reconstruction of the target orbit  corresponding to the initial conditions. From the numerical results of the convergence region analysis on initial state errors, it can be concluded that both $L_1$ and $L_2$ Lyapunov orbits exhibit favorable overall convergence properties, More precisely, $[-0.5\dot{y}_0, 0.21\dot{y}_0]$ for $L_1$ Lyapunov orbit and $[-0.5\dot{y}_0, 0.5\dot{y}_0]$ for $L_2$ Lyapunov orbit. For certain DRO, the convergence domain is relatively limited, $[-0.35\dot{y}_0, 0.5\dot{y}_0]$ for low-amplitude DRO and $[-0.09\dot{y}_0,$ $  0.3\dot{y}_0]$ for medium-amplitude DRO. Furthermore, Halo orbits show a well-defined convergence region $[-0.5\dot{y}_0, 0.4\dot{y}_0]$ across the entire optimizable parameter range. The above numerical findings offer a general reference for other initial state perturbation ranges of the orbit types considered;  however, they may not be strictly identical. Thus, further investigation and analysis are required to determine the influence factors for the convergence domains for all orbital families.

\section*{Funding  }
This work is supported by the National Natural Science Foundation (NNSF) of China under Grant Nos. 12322202, 12172013 and 12371386.
\section*{Author contribution}
Y. Chen: Formal analysis, Methodology, and Writing – original draft; M. Cui: Funding acquisition; Y.J. Qian: Funding acquisition and Supervision; Y.F. Cheng: Investigation; W.X. Chen: Writing – review and editing.


\begin{thebibliography}{[34]}
%
% and use \bibitem to create references. Consult the Instructions
% for authors for reference list style.
%
%\bibitem{RefJ}
% Format for Journal Reference
%Author, Article title, Journal, Volume, page numbers (year)
% Format for books
%\bibitem{RefB}
%Author, Book title, page numbers. Publisher, place (year)
% etc
	\setcounter{enumiv}{0}
	\bibitem{HKC2001}Howell, K.C. 
Families of Orbits in the Vicinity of the Collinear Libration Points[J]. 
J Astronaut Sci. \textbf{49(1)}, 107-125 (2001) %,doi:\href{https://doi.org/10.1007/BF03546339}
	
	
	\bibitem{CELLETTI2025}Celletti, A., {De Blasi}, I.  and {Di Ruzza}, S.
	Perturbative methods and synchronous resonances in Celestial Mechanics.
	Appl. Math. Model. \textbf{143}, 116040 (2025)
	
    %,doi: \href{https://doi.org/10.1007/s10509-017-3115-y}.



	\bibitem{HKC1984}Howell, K.C. 
Three-dimensional, periodic, 'halo' orbits. Celest Mech Dyn Astr. \textbf{32(1)}, 53-71 (1984) %,doi:\href{https://doi.org/10.1007/BF03546339}

    \bibitem{QYY2017}Qian, Y.J., Yang, X.D., Zhai, G.Q. and et al..
	Analytical and numerical construction of vertical periodic orbits about triangular libration points based on polynomial expansion relations among directions.
	Astrophys Space Sci. \textbf{362}, 136 (2017)
    
    \bibitem{FK1973}Farquhar, R. and Kamel, A. 
    Quasi-periodic orbit about the translunar libration point.
	Celest Mech Dyn Astr, \textbf{7}, 458-473(1973)
    %,doi:\href{https://doi.org/10.1007/BF01227511}
	
	
    
	\bibitem{RDL1980}Richardson, D.L.
	Analytic Construction of Periodic Orbits about the Collinear Points.
	Celest Mech Dyn Astr. \textbf{22(3)}, 241–253 (1980)
    %,doi: %\href{https://doi.org/10.1007/BF01229511} 

	\bibitem{PENG2017} Peng, H., X.L Bai and S.j Xu
	Continuation of periodic orbits in the Sun-Mercury elliptic restricted three-body problem.
	Commun nonlinear SCI. \textbf{47}, 1-15 (2017)
    %,doi: %\href{https://doi.org/10.1007/BF01229511} 

    \bibitem{QYZ2018}Qian, Y.J., Yang, X.D., Zhang, W. and Zhai, G.Q.
	Periodic Motion Analysis Around the Libration Points by Polynomial Expansion Method in Planar Circular Restricted Three-body Problem
	Nonlinear Dyn. \textbf{91(1)}, 39-54 (2018)

    \bibitem{QYY2018}Qian, Y.J., Yang, L.Y., Yang, X.D. and et al.
	Parametric stability analysis for planar bicircular restricted four-body problem. 
	Astrodyn. \textbf{2}, 147–159 (2018). 
	%,doi:\href{https://doi.org/10.1007/s42064-017-0017-2}

    \bibitem{Y2022}Yuika, K. and Mitsuru, S.
	Variational existence proof for multiple periodic orbits in the planar circular restricted three-body problem.
	Nonlinearity. \textbf{35(3)}, 1431-1446 (2022)

     \bibitem{C2023}Canales, D., Perera, S.M., Kurttisi, A. and et al.
	A Low-Complexity Algorithm to Determine Trajectories Within the Circular Restricted Three-Body Problem. J Astrnaut Sci \textbf{70(6)}, 46-68 (2023)


    \bibitem{F2024} Tagliaferri,F., Emmanuel, B. and Giacomo, A.
	  Global Optimization for trajectory Design via Invariant Manifolds in the Earth-Moon Circular Restricted Three-Body Problem.
	arXiv preprint: 2405.18916 (2024)


    \bibitem{YGI2025}Yamaguchi, K., Gu, X.B., Inamori, T., and et al.
	Feasibility of Orbital Capture of Near-Earth Asteroids Based on the Planar-Circular Restricted Three-body Problem[J]. 
	Adv Space Res, \textbf{75(4)}, 3806-3819 (2025)


    \bibitem{PLH2025} Peng, L, Liang, Y.Y. and He, X.J.
	Transfers to Earth-Moon Triangular Libration Points by Sun-Perturbed Dynamics. 
	Adv Space Res, \textbf{75(3)}, 2837-2855 (2025)

    \bibitem{YG2020}Li, J., Zhang, H., Liu, Y., and  et al. 
	Utilizing symmetry-enhanced physics-informed neural network to obtain the solution beyond sampling domain for partial differential equations. Nonlinear
Dyn. \textbf{111}, 21861–21876 (2023)
	%, doi: \href{https://doi.org/10.1016/j.jcp.2020.109409}{10.1016/j.jcp.2020.109409}

    \bibitem{SAQLAIN2023} Saqlain,S., Zhu, W., Charalampidis,E.G. and Kevrekidis P. G. 
	Discovering governing equations in discrete systems using PINNs.
	Commun nonlinear SCI. \textbf{126}, 107498 (2023) 

    \bibitem{YL2024}Ma, P.L., Pu, J.C. Peng, W.Q. PINN for solving localized wave solutions and inverse problems involving the Gerdjikov-Ivanov equation. Nonlinear Dyn. \textbf{113}, 26547–26559 (2025)


    \bibitem{WANG2025} Wang,C.Y., Shi, F., Di, Y.J. and Liang,H.
	Multi-task learning enhanced deep neural networks for solving delay differential equations. Commun nonlinear SCI. \textbf{}, 109528 (2025)
    
    \bibitem{WANG2026} Wang, R., Hao, X.H.  and Zhang, L.P., Wu, W.c.  and Hou, S.B. 
	Coaxial permanent magnetic coupling analysis based on physics-informed neural network with periodic fixed boundary conditions.
    Appl. Math. Model \textbf{149}, 116301 (2026)

    \bibitem{AR2020} Rubinsztejn, A., Sood, R.
	Neural network optimal control in astrodynamics: Application to the missed thrust problem.
	Acta Astronaut.  \textbf{176}, 192-203 (2020) %,doi:\href{https://doi.org/10.1016/j.actaastro.2020.05.027}
    
    \bibitem{RPK2019} Raissi, M., Perdikaris, P. and Karniadakis, G.E.
	Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations.
	J. Comput. Phys. \textbf{378}, 686-707 (2019)
    %, doi:  \href{https://doi.org/10.1016/j.jcp.2018.10.045}{10.1016/j.jcp.2018.10.045}
	%,doi:\href{https://iopscience.iop.org/article/10.1088/1361-6544/ac4c2b}.


    \bibitem{LLF1998} Lagaris, I.E., Likas, A. and Fotiadis, D.I.
	Artificial neural networks for solving ordinary and partial differential equations.
	IEEE Trans. \textbf{9(5)}, 987-1000 (1998)
	%, doi:  \href{https://ieeexplore.ieee.org/document/712178}{10.1109/72.712178}

    \bibitem{XDX2024} Zhou, X., Qiao, D. and  Li, X.
	Neural Network-Based Method for Orbit Uncertainty Propagation and Estimation.
	T-AES.\textbf{60(1)}, 1176-1193 (2024)  
	%,doi:\href{https://ieeexplore.ieee.org/document/10316560}{10.1109/TAES.2023.3332566}



	
	\bibitem{GMZ2023} Gong, B., Ma, Y., Zhang, W. and et al.
	Deep-neural-network-based angles-only relative orbit determination for space non-cooperative target.
	Acta Astronaut. \textbf{204}, 552-567 (2023)  
	%,doi:\href{https://doi.org/10.1016/j.actaastro.2022.09.024}{10.1016/j.actaastro.2022.09.024}
    
	\bibitem{ZL2024} Zhou, X., Li, X. and Zhang, Z.
	Neural Network–Assisted Initial Orbit Determination Method for Libration Point Orbits. 
	J. Aerosp. Eng. \textbf{37(4)}, 04024046 (2024)  
	%, doi:\href{https://doi.org/10.1061/JAEEEZ.ASENG-5482}{10.1061/JAEEEZ.ASENG-5482}

	\bibitem{EDD2022} Schiassi, E., D’Ambrosio, A., K, Drozd. and et al.
	Physics-informed neural networks for optimal planar orbit transfers.
	J. Spacecr. Rockets. \textbf{59(3)}, 834-849 (2022)
    %,doi:  \href{https://doi.org/10.2514/1.A35138}{10.2514/1.A35138}

    \bibitem{DAE2021}D'Ambrosio, A., Schiassi, E., Curti, F. and Furfaro, R.
	Physics-Informed Neural Networks Applied to a Series of Constrained Space Guidance Problems.
	31st AAS/AIAA Space Flight Mechanics Meeting tenutosi a Virtual. (2021)


	\bibitem{AAS2021} Krishnapriyan, A., Gholami, A. and et al.
	Characterizing Possible Failure Modes in Physics-Informed Neural Networks, Advances[J].
	NeurIPS \textbf{34},26548–26560 (2021) 
	%,doi: \href{https://arxiv.org/abs/2109.01050}{abs/2109.01050}

	\bibitem{KA2022} Haitsiukevich, K. and Ilin, A. 
	A Limited Memory Algorithm for Bound Constrained Optimization.
	arXiv:2204.05108 (2022)
    %,doi:  \href{https://arxiv.org/abs/2204.05108}{abs/2204.05108}

    \bibitem{GYW2023}Guo, J., Yao, Y. and et al.
	Pre-training strategy for solving evolution equations based on physics-informed neural networks[J].
	J. Comput. Phys. \textbf{489}, 112258 (2023)
    %,doi:  \href{https://doi.org/10.1016/j.jcp.2023.112258}{10.1016/j.jcp.2023.112258}


	\bibitem{LP2021} Lu, L.,  Pestourie, R. and et al.
	Physics-Informed Neural Networks with Hard Constraints for Inverse Design.
	SIAM. \textbf{43}, 1064-8275(2021)
    %,doi:\href{https://doi.org/10.1137/21M1397908}{10.1137/21M1397908}.

    
	\bibitem{CBE2025}Bihlo, E.C. and Bihlo, A.
	Exactly conservative physics-informed neural networks and deep operator networks for dynamical systems.
	ANN.\textbf{181}, 106826 (2025)
    %,doi: \href{https://doi.org/10.1016/j.neunet.2024.106826}{10.1016/j.neunet.2024.106826}. 



	\bibitem{RPJC1995}Byrd, R.H.,  Lu, P., Nocedal, J. and Zhu, C.
	Improved training of physics-informed neural networks with model ensembles.
	SIAM J. Sci. Comput. \textbf{16(5)}, 1190–1208 (1995)  
	%.doi:\href{https://doi.org/10.1137/0916069}{10.1137/0916069}


	\bibitem{JMLR2015}Baydin,A.G., Pearlmutter, B.A.,  Radul, A. and Siskind, J.M. 
	Automatic differentiation in machine learning: a survey.
	J. Mach. Learn. Res.. \textbf{18}, 153:1-153:43 (2015)  
	%.doi:\href{https://doi.org/10.1137/0916069}{10.1137/0916069}


	

	



\end{thebibliography}

\end{document}



