Fast and accurate adaptive collocation iteration method for orbit dynamic problems

Honghua DAI, Zhe ZHANG,*, Xuechuan WANG, Haoyang FENG,Changtao WANG, Xiaokui YUE

a National Key Laboratory of Aerospace Flight Dynamics (AFDL), Northwestern Polytechnical University, Xi’an 710072, China

b School of Astronautics, Northwestern Polytechnical University, Xi’an 710072, China

KEYWORDS

Abstract For over half a century,numerical integration methods based on finite difference,such as the Runge-Kutta method and the Euler method, have been popular and widely used for solving orbit dynamic problems.In general, a small integration step size is always required to suppress the increase of the accumulated computation error, which leads to a relatively slow computation speed.Recently, a collocation iteration method, approximating the solutions of orbit dynamic problems iteratively, has been developed.This method achieves high computation accuracy with extremely large step size.Although efficient,the collocation iteration method suffers from two limitations: (A) the computational error limit of the approximate solution is not clear; (B) extensive trials and errors are always required in tuning parameters.To overcome these problems, the influence mechanism of how the dynamic problems and parameters affect the error limit of the collocation iteration method is explored.On this basis, a parameter adjustment method known as the‘‘polishing method” is proposed to improve the computation speed.The method proposed is demonstrated in three typical orbit dynamic problems in aerospace engineering: a low Earth orbit propagation problem,a Molniya orbit propagation problem,and a geostationary orbit propagation problem.Numerical simulations show that the proposed polishing method is faster and more accurate than the finite-difference-based method and the most advanced collocation iteration method.

1.Introduction

Orbit dynamic equations are important in aerospace engineering.1–5Since the majority of these equations cannot be solved analytically, numerical integration methods are essential in practical problems.Traditionally, Orbit dynamic problems are always solved by finite-difference-based methods, such as the Euler method and the Runge-Kutta method.However,these methods demand very small step sizes to achieve high accuracy, which makes them inefficient.In order to overcome this issue, researchers developed a series of methods of better performance.

In the last two decades, a collocation iteration method has become a hot topic.This method solves dynamic problems by iteratively revising approximate solutions,and shows excellent efficiency as well as accuracy.6–17The idea of solving dynamic equations through iteration originates from Picard’s research,wherein the Picard method was developed by combining the iteration method and the definite integral formula.18Though the iteration formula is concise,the Picard method is not practical in engineering since it requires extremely complex symbolic calculations.Clenshaw and Norton19combined the Picard method with the collocation method,and proposed the first collocation iteration method, which is radically different from finitedifference-based methods.This modification frees the Picard method from tedious symbolic calculations,but is not accurate and there is a complex calculation process in each iteration.To simplify and accelerate Clenshaw’s method, the Chebyshev-Picard method using the Chebyshev polynomial approximation to obtain the global solutions of dynamic problems was proposed by Fukushima.6He introduced this method into aerospace engineering and boosted the subsequent studies.Inspired by Fukushima’s research, a vector–matrix form of the Chebyshev Picard method named the Modified Chebyshev-Picard Iteration method (MCPI) was proposed by Bai and Junkins.7Bai8solved a series of orbit dynamic problems using the MCPI,and her research showed that this method is significantly faster and more accurate than traditional finitedifference-based methods.To further enhance the performance of the MCPI,Wang et al.9presented the Local Variational Iteration Method (LVIM) by combining the variational iteration method and the collocation method.Their research showed that the MCPI is a special case of the LVIM,which holds advantages over the MCPI in terms of efficiency and accuracy in orbit dynamic problems.The aforementioned methods are efficient tools in solving ordinary differential equations,but cannot handle two-point boundary value problems directly.To further expand the application scope of the collocation iteration method, Feng et al.20proposed the Quasi-linear Local Variational Iteration Method (QLVIM).Numerical simulations showed that the QLVIM can solve two-point boundary value problems with high efficiency.All the aforementioned methods constitute a complete system of high performance tools for solving nonlinear orbit dynamic problems.

Though the collocation iteration method has the potential to compute with large step sizes, its efficiency is sensitive to system parameters such as the step size and the number of collocation points.This shortcoming hinders the dissemination of this kind of method,since it always requires extensive trial and error to find proper parameters before computation.To solve this problem,Woollands and Junkins16presented the Adaptive Picard-Chebyshev numerical integration method (APC).This method divides the time domain of the perturbed two-body problem into several segments with the same true anomaly,and increases the number of collocation points in each segment to satisfy the user-specified precision tolerance.Numerical simulations showed that this method can successfully decrease the computation time.However, the APC requires prior knowledge of the solution before computation;thus,it is not suitable for the problems with little prior knowledge.To overcome this limitation,Wang et al.21proposed another self-tuning collocation iteration method, ph-LVIM, which adjusts parameters according to the computation error22–28and achieves high performance in solving orbit propagation problems.Unfortunately, it may get stuck in an adaptive iteration.The two parameter-adjustment methods provide different ideas to enhance the performance of collocation iteration methods,but both have significant drawbacks.

To further utilize the potential of collocation iteration methods, we explore how parameters affect the convergence speed of collocation iteration methods through theoretical analysis and numerical illustrations.After that, we present an adaptive parameter adjustment method which divides intervals at the places where collocation iteration methods converge slowly.Our method is validated in three benchmark problems in astrodynamics: a Low Earth Orbit (LEO) propagation problem,a Molniya orbit propagation problem,and a Geostationary Orbit (GEO) propagation problem.

The main contributions of this paper are threefold.First,the influence of the second order derivative and the integration matrix on the convergence speed of collocation iteration methods is revealed.Second, an adaptive collocation iteration method which can greatly improve the computational performance of existing methods is presented to solve the nonlinear dynamic problems efficiently and precisely.Third,a new strategy to set the initial value by combining cold start and warm start is presented to further increase the efficiency without adding extra computation load.

This paper is organized as follows.In Section 2, some typical collocation iteration methods are reviewed.In Section 3,the convergence speed of typical collocation iteration methods is analyzed.The analysis in Section 3 is validated by numerical simulations in Section 4.In Section 5, the parameter adjustment method is presented.In Section 6, the performance of the proposed method is compared with two other adaptive methods in three typical orbit propagation problems.This work is concluded in Section 7.

2.Typical collocation iteration methods

The collocation iteration method is a numerical method for solving differential equations,which is radically different from the traditional step-by-step finite-difference-based methods.In this section, we introduce two typical collocation iteration methods that will be used in the rest of this paper.

2.1.MCPI method

The modified Chebyshev-Picard iteration method is a collocation iteration method of vector–matrix form combining the Picard method,the Chebyshev approximation and the collocation method.

Considering a nonlinear ordinary differential equation with an initial condition

The Picard method solves it through the following iterative integration formula

In practice, the integral operation of the Picard method becomes extremely complex after the first few iterations.In order to simplify the iteration process of the Picard method,the MCPI approximates the x and the ˙x in Eq.(1) by two Nth order orthonormal polynomials

where Tiis the orthonormal polynomial of degree i,and βiand Fiare the coefficients of Ti.This modification transforms the update of the approximate solution in each iteration into the update of the coefficients of the orthonormal polynomials.

The coefficients in Eqs.(3)and(4)are calculated by the following formulas in each iteration8

Substituting Eqs.(3)–(8) into Eq.(2), the iteration formula of the MCPI can be expressed as a vector–matrix form8

2.2.Local variational iteration method

The local variational iteration method combines the variational iteration method and the collocation iteration method.Different from the MCPI, the LVIM uses a feedbackweighted optimal error to correct the approximate solution in each iteration.29

The iteration formula of the variational iteration method is

where λ(τ) is the generalized Lagrange multiplier.

Substituting Eqs.(3)and(4)into Eq.(10)and approximating the generalized Lagrange multiplier by a truncated polynomial series,29the iteration formula of the LVIM can be expressed in the vector–matrix form

where T0is the 0th order differential transformation of the generalized Lagrange multiplier, P is the integration matrix,and D is the differential matrix.

3.Convergence speed of collocation iteration methods

In this section, the convergence speed of the two collocation iteration methods mentioned in Section 2 is analysed theoretically.For convenience, we first analyse the convergence speed of the Picard method and the variational iteration method in a continuous system.After that, we analyse the convergence speed of the MCPI and the LVIM.

3.1.Convergence speed of integration iteration methods in continuous systems

3.1.1.Convergence speed of Picard method

According to Eqs.(1) and (2), the iteration formula of the Picard method satisfies the following equation

Subtracting Eq.(12) from Eq.(2), the following inequality can be derived according to the absolute value inequality

For the orbit propagation problems defined in the Earthcentered Earth-fixed coordinate system, it can be seen from the dynamic equations that the second order derivatives of the state variables exist, and are bounded at any point except the origin.Therefore, there exists a constant M satisfying the following formula

According to Eq.(14),for any two points(t,x1)and(t,x2)in the rectangular domain[t0,tf]×[x0-δ,x0+δ],function f satisfies the Lipschitz condition

where L is a Lipschitz constant which is no greater than M.

Denoting the difference between the nth and the (n - 1)th approximations of x as dn

substituting Eq.(12) into Eq.(16), we have

Denoting the maximum value of dnin the domain[t0,tf]as

According to Eqs.(13), (16) and (18), we have

According to Eq.(19),dncan be estimated by dn-1.Repeating the amplification of di(i = n, n - 1,???,1) in Eq.(19), we have

A necessary condition for xnto be convergent is

According to Eq.(20),a sufficient condition for Eq.(21)to hold is

that is

According to Eq.(16), the nth approximation of x can be represented by di(i = 1,2,???,n) and the initial guess x0

Assuming the exact solution for Eq.(1) is x*, the upper bound of the error between x*and xncan be derived according to the formula for the sum of a geometric series

Without loss of generality, we can assume t = 1 and t0= -1.The upper bound of the computation error given by Eq.(25) is shown in Fig.1.It can be seen that the smaller the Lipschitz constant is,the fewer the iterations are needed to reach the given tolerance.Though the Lipschitz constants in many problems are not known, they are no bigger than the maximum value of the second derivative of the exact solution according to Eqs.(14) and (15).Approximating the Lipschitz constant in Eq.(25) by the maximum value of the second derivative of the solution, we can draw the conclusion that the Picard method converges fast where the upper bound of the second derivative of the solution is small.

3.1.2.Convergence speed of variational iteration method

Replacing the generalized Lagrange multiplier in Eq.(10) by its 0th order Taylor series approximation, the iteration formula of the variational iteration method can be expressed as follows:29

According to Eq.(26), we have

Fig.1 Upper bound of approximation error with different Lipschitz constant.

Following the same analysis procedure as that of the Picard method, we can draw the conclusion that the variational iteration method converges fast where the upper bound of the second derivative of the solution is small.

3.2.Convergence speed of integration iteration methods in discrete systems

The numerical solutions obtained by collocation iteration methods are discrete; therefore, the analysis procedure in Section 3.1 does not fully apply to collocation iteration methods.In this part,the convergence speed of the two collocation iteration methods mentioned in Section 2 is analysed specifically.

3.2.1.Convergence speed of modified Chebyshev Picard iteration method

Since the integrations of Chebyshev polynomials can be calculated by their properties,9the iteration formula of the MCPI can be written as the following simplified form

where P is the integration matrix.

We can obtain the following formulas according to Eqs.(15) and (28)

Denoting the difference between the nth and the (n - 1)th approximation of x as dn, the upper bound of the approximation error can be derived as

3.2.2.Convergence speed of local variational iteration method

SimilartotheanalysisprocedureoftheMCPIinSection3.2.1,we can obtain the following formula according to Eqs.(11)and(15)

It can be seen from Eqs.(30)and(32)that the convergence speed of collocation iteration methods is closely related to the Lipschitz constant and the F-norm of the integration matrix.Similar to the iteration methods analyzed in Section 3.1, we can draw the conclusion that the smaller the upper bound of the second derivative of the exact solutions and the F-norm of the integration matrices are, the higher the convergence speed of the two collocation iteration methods is.

3.3.F-norm of integration matrix

Section 3.2 shows that the F-norm of the integration matrix affects the convergence speed of collocation iteration methods.In this section,we search for the best approximate polynomials according to the F-norm of the integration matrix.

Legendre polynomials, Chebyshev polynomials and trigonometric polynomials are defined in finite intervals, such as[-1,1]or[-π,π],so they can be used in collocation iteration methods to solve problems in finite intervals.The F-norms of the integral matrices of these three polynomials of different degree are shown in Fig.2.The collocation points of the integration matrices are selected as Chebyshev-Gauss-Lobatto(CGL) nodes.

In Fig.2,we can see that the F-norms of the integral matrices of Chebyshev polynomials are always the smallest among the three orthogonal polynomials.According to Eqs.(30)and(32),Chebyshev polynomials are the best orthogonal polynomials to improve the convergence speed of collocation iteration methods.

4.Verification of theoretical analysis

In this section,the conclusion obtained in Section 3 that collocation iteration methods converge fast in intervals with small second derivative of the solution is verified in two numerical simulations.

Fig.2 F-norm of the integral matrices with different number of collocation points.

4.1.Verification through ph-LVIM

The ph-LVIM adjusts parameters according to the computation error and the convergence speed.If the algorithm converges too slow or the error exceeds a prescribed tolerance,the ph-LVIM will shorten the step size and increase the number of collocation points in the subinterval to ensure that the prescribed tolerance can be attained.Since the ph-LVIM can adaptively decrease the step size and increase the number of collocation points in the subintervals where the approximate solution cannot attain the prescribed accuracy within a given time, we use it to monitor the convergence speed of the collocation iteration method.

To verify the conclusion obtained in Section 3, we use ph-LVIM to solve the orbit propagation problem that the norm of the second derivative of the solution changes significantly.It is solved by the ph-LVIM with two sets of initial parameters.Through many trials, we obtain two sets of initial parameters that enable ph-LVIM to adjust the number of collocation points without changing the step size in the calculation process.These simulations will clearly show the domain which is difficult to converge by reducing variables.Solutions and parameters of the ph-LVIM are shown in Figs.3-6.The ordinates in Figs.4 and 6 show the number of collocation points in each subinterval.Curves representing norms of velocity(the first derivative)and acceleration(the second derivative)in Figs.4 and 6 have been compressed and panned to fit the figures.

As shown in Figs.4 and 6,the step sizes in the two simulations keep constant, and the number of the collocation points is adjusted automatically in each subinterval.It is worth noting that the trend of the number of the collocation points and the second derivative of the solution are basically the same.This fact is in good agreement with the conclusion we draw in the last section that the collocation iteration method requires more iterations to achieve the prescribed accuracy in the subintervals with large second derivatives.

However,we can also see that the trends of the first and the second derivatives of the solution for this problem are exactly the same.To exclude the interference introduced by the first derivative, we further carry out the following simulation.

Fig.3 Orbit obtained via Parameter set 1.

Fig.4 Step size, velocity, acceleration of solution and approximate order of ph-LVIM.

Fig.5 Orbit obtained via Parameter set 2.

Fig.6 Step size, velocity, acceleration of solution and approximate degree of ph-LVIM.

4.2.Verification through collocation iteration method with fixed parameters

In this section, several orbits that have inverse relationships between the magnitude of their first and second derivatives are solved by a collocation iteration method with constant parameters.The convergence speed of the collocation iteration method in different subintervals are compared through the number of iterations required to attain the prescribed accuracy.

The solutions for these orbit propagation problems have the following characteristics:

(1) The 2-norm of the first derivative of the orange orbit is smaller than that of the black orbit near the perigee.

(2) The 2-norm of the first derivative of the green orbit near the perigee is greater than the 2-norm of the first derivative of the blue orbit.

(3) The 2-norm of the first derivative of the magenta orbit is greater than that of the three elliptical orbits near their apogees.

All the elliptical orbits are integrated from their perigees,and the initial values of different orbits are shown in Table 1.The solutions for these orbit propagation problems obtained by the collocation iteration method are shown in Fig.7.The number of iterations required to achieve the prescribed accuracy in different subintervals are shown in Figs.8(a)–(f).

It can be seen in Figs.8(b), (d) and (e) that the number of the iterations required to achieve the prescribed accuracy in a subinterval decreases, when the subinterval gradually moves away from the perigee of the orbit.This means that the convergence speed of the collocation iteration method will increase,when the 2-norm of the first and the second derivatives decrease.Figs.8(a)and 8(d)show that the orange orbit always requires more iterations to achieve the prescribed accuracy than that required in solving the black orbit.Since the 2-norm of the first derivative of the solution near the perigee of the black orbit is greater than that of the orange orbit and the 2-norm of the second derivative of the solution of the orange orbit is always greater, we come to the conclusion that it is the second derivative, instead of the first derivative,that affects the convergence speed of the collocation iteration method.Likewise, though the 2-norm of the first derivative of the solution at the apogee of the magenta orbit is greater than those of the green, black, and red orbits, the 2-norm of the second derivative of the magenta orbit is the smallest among them, and always takes more iterations for the three elliptical orbits than the magenta orbit to meet the prescribed accuracy.These facts further validate the conclusion that the convergence speed of the collocation iteration method will be high in those intervals with small second derivative of the solution.

5.Polishing method

Inspired by the conclusion that the second derivative of the solution has an important impact on the efficiency of collocation iteration methods, we design the following adaptive parameter adjustment method.

As shown in Fig.9, CGL collocation points are sparse at the middle,and dense at the boundary of a domain.Since adding collocation points in a subinterval can reduce the computation error27, we divide the interval at the position where the local maximum point of the second derivative of the approximate solution occurs, and add collocation points nearby.

Table 1 Initial values of orbit propagation problems.

Fig.7 Orbits with different initial states.

Specifically, the time domain is kept intact, and an upper limit of the number of the collocation points is provided first.The number of the collocation points in the domain is set to the upper limit.Then, the differential equations are solved by the collocation iteration method with a cold start, and the computation error will be evaluated after each iteration.If the prescribed accuracy is achieved and the computation error rises no more than once in a row, the result will be recorded and the iterative procedure will be terminated.Otherwise, if the computation error rises twice successively, parameters will be adjusted by the following procedure.

When the computation error rises more than a given limit successively, the approximate solution and its largest second derivative will be recorded.Then, parameters will be adjusted to achieve the prescribed accuracy.If the step size equals to or is smaller than a given lower limit of the interval length, the parameters will be adjusted according to the computation error:

(1) When the method diverges since the first iteration, the number of the collocation points in the interval will be increased by 10%,and the approximate solution at each collocation point will be recalculated with a cold start.That is to say,the initial approximation at each new collocation point is the same as the initial value.

(2) When the method converges in the first few iterations,the number of the collocation points in the interval will be increased by 10%, and the approximate solution at each collocation point will be recalculated with a warm start.That is to say, the initial approximation at each new collocation point is approximated by interpolating the approximate solution obtained in the last iteration.

If the step size is larger than the given lower limit of the interval length, parameters will be adjusted according to both the step size and the number of local maximum point of the second derivative:

(1) If the method diverges since the first iteration,the interval will be divided into two subintervals at the midpoint.The number of the collocation points in each subinterval will be halved,and the approximate value at each collocation point in the new subintervals will be approximated iteratively with a cold start.

(2) If the method converges in the first few iterations but no local maximum point of the second derivative of the approximate solution appears in the interval, the interval will be divided into two subintervals at the midpoint.The number of the collocation points in each subinterval will be halved, and the solution at each collocation points in the new subintervals will be approximated iteratively with a warm start.

(3) If there is only one local maximum point of the second derivative of the solution in the interval,the interval will be divided at the position where the local maximum point occurs.Then, the value at each collocation point in the subinterval will be recalculated with a warm start.

(4) If there is more than one local maximum point of the second derivative of the solution in the interval, the interval will be divided into two subintervals at the positions where the largest and second largest local maximum point occur.Then, the value at each collocation point in the subinterval will be recalculated with a warm start.

In all the four cases mentioned above, the number of the collocation points in each subinterval is proportional to the interval length.If the number of the collocation points is less than the given lower limit to the number of the collocation points in any subinterval after division, it will be increased to the lower limit immediately.

The preceding mentioned parameter adjustment process will be repeated until the approximate solutions in all intervals satisfy the prescribed precision tolerance.

Our method always divides the intervals at the position of the maximum second derivative of the solution,which reminds us that carpenters always begin to polish a board from the roughest place, and then move to the less rough places until the material is as smooth as required.Therefore, we name it the ‘polishing method’.

Fig.8 Iteration number of different orbits.

6.Numerical simulations

In this section, the effectiveness of the polishing method combining with the second order LVIM29is validated in solving three typical orbit dynamic problems: a Low Earth Orbit(LEO) propagation problem, a Molniya orbit propagation problem,and a Geostationary Orbit(GEO)propagation problem.The numerical simulation is carried out in MATLAB R2020b, using a computer with Intel(R) Core(TM) i7-10700F 2.9 GHz CPU, 32 G RAM.The performance of the polishing method is compared with that of ode45, the implementation of DP5(4) in MATLAB, and ph-LVIM, the most advanced adaptive collocation iteration method.Parameters,initial values and the computation time of different problems are shown in Tables 2-10.The orbits obtained by different methods are shown in Figs.10-12.The benchmarks in all the three simulations are obtained by setting the Rel.Error and Abs.Error of ode45 to 10-16and 10-20,below the lowest value allowed by MATLAB.Computation errors of different methods are shown in Figs.13-15.

Fig.9 Distribution of CGL points.

Table 2 Initial value of LEO.

Table 3 Computation time of different methods for LEO.

Table 4 Parameters of different methods for LEO.

Table 6 Computation time of different methods to solve Molniya orbit.

Table 7 Parameters of different methods for Molniya orbit.

Table 8 Initial value of GEO.

Table 9 Computation time of different methods for GEO.

Table 10 Parameters of different methods for GEO.

Fig.10 LEO obtained via different methods.

Fig.11 Molniya orbit obtained via different methods.

Fig.12 GEO obtained via different methods.

In Fig.13, the computation error of the solution for the LEO propagation problem obtained by the polishing method,the ph-LVIM and ode45 are compared.It can be seen that the computation error of all the three methods increases fast at the beginning.The error of the polishing method is roughly the same as that of the ph-LVIM in the first 1000 s.After that,the error of the ph-LVIM gradually overtakes that of the polishing method, and the computation error of the polishingmethod becomes the smallest among the three methods.It can be clearly seen that the computation error of the other two methods increases faster than that of the polishing method.This is because the polishing method can preferentially select larger step size to reduce the accumulated error caused by too many subintervals.At the end of the time domain, the error of the polishing method is more than one order of magnitude smaller than that of the other two methods.As larger step size is always used in the polishing method,the computation load will be significantly reduced, and the computation efficiency will be improved.Table 3 shows that the polishing method is 140% faster than the ph-LVIM, and 305% faster than ode45 in solving the LEO propagation problem.

Fig.13 Position error of LEO obtained via different methods.

Fig.14 Position error of Molniya orbit obtained via different methods.

Fig.15 Position error of GEO obtained via different methods.

In Fig.14, the computation error of the solution for the Molniya orbit propagation problem obtained by the three methods are compared.It can be seen that the computation error of the solution obtained by all the three methods basically increases with time in this simulation.Though the computation error of the polishing method overtakes that of the ode45 for about 5000 s after 22000 s, it is the smallest one among those of the three methods during most of the time domain.For the same reason as the last simulation that the polishing method will significantly decrease the computation load, Table 7 shows that the polishing method is 156% faster than the ph-LVIM, and 28% faster than ode45 in solving the Molniya orbit propagation problem.

In Fig.15, the computation error of the solution for the GEO propagation problem obtained by the three methods are compared.Due to less computation load and the computation error caused by larger step size, it can be seen that the computation error of the polishing method has the slowest growth rate, and is the smallest one among the three methods for most of the time.Table 9 shows that the polishing method is 311% faster than the ph-LVIM, and 91% faster than ode45 in solving the GEO propagation problem.

In sum,the polishing method is the most accurate and efficient method compared with ode45 and the ph-LVIM in the three simulations.

7.Conclusions

This paper develops an adaptive parameter adjustment method for collocation iteration methods.The convergence speed of typical collocation iteration methods is analyzed through theoretical analysis and numerical simulations.The second derivatives of solutions and the norms of integration matrices are found to significantly affect the convergence speed of collocation iteration methods.Based on these findings, the polishing method that can adaptively divide intervals at the places where the local maximum point of the second derivative of the solution occurs is proposed to accelerate the convergence speed of collocation iteration methods.The effectiveness of the proposed method is verified in solving three orbit dynamic problems.Numerical simulations demonstrate that the polishing method is at least 28% faster than ode45, and 140% faster than the ph-LVIM.It is also the most accurate method among them.The proposed method is highly efficient,and can be used in a series of nonlinear dynamic problems where accurate solutions are sought within a very limited time.

Finally, we point out that though the polishing method is very efficient and the most precise method when compared with traditional methods, it may loss some precision at the beginning since it always uses more sparse collocation points than the other methods.In the future, we will further explore this problem and try to overcome it.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgements

This study was co-supported by the National Key Research and Development Program of China (No.2021YFA0717100)and the National Natural Science Foundation of China(Nos.12072270 , U2013206).