Yang Hong,Li Xin-hong,Ding Wen-zhev

Department of Aerospace Science and Technology,Space Engineering University,No.1 BaYi Road,HuaiRou District,Beijing,101416,China

Keywords:Tangent impulse interception Minimum time Interception uncertainty Multi-objective optimization Hybrid optimization

ABSTRACT The traditional tangent impulse interception problem does not consider the in fluence of actual deviation.However,by taking the actual state deviation of the interceptor into the orbit design process,an interception orbit that is more robust than the nominal orbit can be obtained.Therefore,we study the minimum time interception problem and the minimum terminal interception error problem under tangent impulse conditions and give an orbit optimization method that considers the interception time and the interception uncertainty.First,we express the interceptor's transfer time equation as a form of flight path angle,establish a global optimization model for solving the minimum time tangent impulse interception and give a hybrid optimization algorithm based on Augmented Lagrange Genetic Algorithm-Sequential Quadratic Programming(ALGA-SQP).Secondly,we use the universal time equation and Bootstrap resampling technology to calculate the interceptor's terminal error distribution and establish the relevant global optimization model by using the circumscribed cuboid volume of the interceptor's terminal position error ellipsoid as the optimization index.Finally,we combined the above two singleobjective optimization models to establish a global multi-objective optimization model that considers interception time and interception uncertainty and gave a hybrid multi-objective optimization algorithm based on Non-dominated Sorting Genetic Algorithm II-Goal Achievement Method(NSGA2-GAM).The simulation example veri fies the effectiveness of this method.

1.Introduction

With the development of light gas propulsion technology,the traditional large-scale ground light gas energizing device has successfully achieved miniaturization[1].The device can be used to launch small unguided interceptors at high speed,which provides a new solution for space interception tasks such as debris removal.Since the tangent interception impulse is simple in impulse direction,it is bene ficial for the spacecraft to adjust the attitude at the impulse point.Therefore,the tangent impulse is suitable for use in space with the above device.In view of the speci fic advantages of the tangent impulse,there have been many related studies on it in recent years.Adamyan et al.studied the bitangent transfer problem and gave the corresponding analytical solution through the geometric method[2].Quarta et al.studied how to use tangent impulses for optimal transfer between coplanar ellipses,and gave a simple solution[3].However,the above research cannot be used for space impulse interception problems due to the lack of flight time restrictions on the interceptor and the intercepting target.For the tangent impulse interception and tangent impulse rendezvous that require the interceptor to have the same flight time as the target,Zhang et al.derived the existence conditions of the interception solutions under three transfer methods:tangent to the initial orbit,tangent to the target orbit,and bitangent to the initial orbit and the target orbit[4].At the same time,they also studied a series of tangent impulse problems such as tangent rendezvous between elliptical orbits[5]and tangent rendezvous between the elliptical orbit and hyperbolic orbit[6].

To avoid unexpected situations during the interception process,it is necessary to implement a shorter time interception on the basis of tangent interception.For this,Zhang et al.used iterative calculation methods to carry out relevant research on the problem of minimum time interception under tangent impulse conditions[7].However,the iterative calculation method adopted not only consumes a lot of time and computing resources but also is relatively cumbersome in actual execution.In this regard,Ding et al.used optimization methods to study the minimum time tangent impulse interception problem under relative motion conditions[8].Compared with the iterative calculation method,the optimization method can search for a better solution that meets the accuracy requirements more quickly and conveniently.According to the different search scopes,optimization methods can be divided into global optimization and local optimization.The global optimization algorithm mainly includes particle swarm optimization[9],ant colony algorithm[10],genetic algorithm[11],gray wolf algorithm[12],whale algorithm[13],fruit fly algorithm[14],etc.Its advantage is that it can obtain a relatively optimal solution in the global scope without providing the initial search value.Traditional local optimization algorithms such as the interior point method[15],sequence quadratic programming method [16],trust region method[17],etc.need to provide initial search values,but they can search for local optimal solutions more quickly and accurately.Therefore,we combine the advantages of the above two methods,and use the search results of the global optimization algorithm as the initial search value of the local optimization algorithm,so as to quickly search for a global minimum time tangent impulse interception solution comparable to the iterative calculation results.Considering that the minimum time tangent impulse interception problem includes non-linear constraints such as integer constraints and inequality constraints,we use the genetic algorithm that incorporates the augmented Lagrangian method as the global optimization algorithm.At the same time,to ensure a high local search speed,we adopt the SQP algorithm based on gradient calculation as the local optimization algorithm.

The above research is mainly designed for the nominal orbit,but the actual interception process will be affected by a series of errors such as orbit determination error and impulse error.This will cause the real orbit of the interceptor to deviate from the nominal orbit,which will affect the interception result of the interceptor at the terminal position.Therefore,if the actual state deviation of the interceptor is taken into account in the design process of the orbit,an interceptor orbit that is more robust than the nominal orbit can be obtained.For this,Luo et al.used the linear covariance propagation method(LINCOV)to study the robust optimization problem of multi-impulse rendezvous under the condition of fixed initial impulse point[18].And the in fluence of impulse number and impulse speed on the final rendezvous error of spacecraft is analyzed.On this basis,Yang et al.studied the robust optimization problem of multi-impulse rendezvous considering the J,J,Jperturbation,and orbit re-planning process under the condition of fixed initial impulse point based on the unscented transformation(UT)method[19].In the above studies,the rendezvous orbit types are all set to ellipse,but due to different impulse conditions,the shape of the rendezvous orbit may be hyperbolic or parabolic in addition to the ellipse[20].To reduce the calculation dif ficulty of deviation propagation,we established the universal time equations for unifying three different conic orbits.In this way,by calculating the universal variable at the interception point,we can directly obtain the state of the interceptor at the interception point from the Lagrangian coef ficient.And to make the calculation results more accurate,we use the Bootstrap resampling technology to calculate the terminal error distribution of the interceptor.Besides,since small unguided interceptors cannot maneuver during the interception process,the entire interception process can only be determined by the initial impulse state,so we cannot reduce the terminal error by optimizing the impulse velocity as in the above studies.To this end,we nested the deviation calculation method into the hybrid optimization algorithm ALGA-SQP and optimized the position of the initial impulse point to reduce the terminal position error of the interceptor.

The above two aspects of the research have optimized the interception time and interception uncertainty.The optimization of interception time can make the interceptor hit the target quickly,and the optimization of robustness can make the interception error of the interceptor when it reaches the predetermined interception point smaller.However,in the whole interception process,the total interception time of the interceptor and the terminal interception error are naturally intrinsically related.Longer interception time will lead to serious divergence of terminal errors when the interceptor reaches the predetermined interception point.The larger initial impulse required for a shorter interception time will bring a larger initial error,which may instead make the terminal error of the interceptor more divergent when it reaches the predetermined interception point.Therefore,it is more practical to integrate interception time and interception uncertainty for uni fied optimization.Compared with other mathematical programming methods,evolutionary algorithms are more suitable for solving multiobjective optimization problems.Among them,NSGA2 uses the fast non-dominated sorting method to obtain the optimal Pareto solution set through crowding degree calculation and elite retention strategy[21].It is widely used in many scienti fic and engineering fields due to its advantages of simplicity and high ef ficiency[22].To further improve the algorithm performance,a large number of scholars have improved the NSGA2 in recent years.Deb et al.proposed an NSGA2 algorithm based on reference points to improve the high-dimensional optimization capability of the algorithm[23].Shim et al.combined non-dominated sorting with target decomposition,thereby improving the optimization performance of the algorithm[24].Qiu et al.proposed an adaptive crossdifferential evolution operator for multi-objective optimization[25].These studies have improved NSGA2 in some respects,but due to the inherent defects of the global optimization algorithm,NSGA2 still lacks in local search capabilities and algorithm convergence.So we propose a hybrid multi-objective optimization algorithm NSGA2-GAM and use the algorithm to solve the problem.The hybrid algorithm adds the optimal solutions of two single objective functions,minimum interception time and minimum terminal position error,to the initial population,thereby expanding the search space of the NSGA2 algorithm and enhancing the global search capability of the algorithm.At the same time,according to the searched Pareto frontier,the hybrid algorithm also uses the goal achievement method(GAM)to convert the multi-objective optimization function into a hybrid function.In this way,the SQP algorithm is further used to solve the converted problem,making up for the lack of the local search ability of the NSGA2 algorithm.

In summary,based on the background of the space tangent impulse interception task,this paper studies the minimum time interception problem and the minimum terminal interception error problem.By combining the above two single-target optimization tasks,a multi-target orbit optimization method that takes into account the intercept time and the intercept uncertainty is proposed.The main contributions of this paper are as follows:

(1)A hybrid optimization algorithm ALGA-SQP is given to solve the problem of minimum time tangent impulse interception.This algorithm combines the advantages of ALGA and SQP,and can quickly search for a more accurate global optimal solution without providing initial values.Compared with the iterative calculation method,using ALGA-SQP to solve the minimum time tangent impulse interception problem can search for a globally optimal solution with the same accuracy more quickly and conveniently.

(2)Based on the universal time equation and Bootstrap resampling,a convenient calculation method that can be applied to solve the propagation of deviation under different shape orbits is given.And by nesting the deviation calculation method into the hybrid optimization algorithm ALGA-SQP,the optimization problem of the minimum terminal interception error of the interceptor is solved.The universal time equation uses universal variables to unify the three different conical orbits,solves the divergence problem of the initial value calculation of the hyperbolic orbit when the eccentricity is small and reduces the dif ficulty of calculating the terminal error of the interceptor.The bootstrap resampling can make the statistical results of terminal errors more stable,thereby reducing Monte Carlo sampling times and saving calculation costs.And by nesting the deviation calculation method into the hybrid optimization algorithm ALGA-SQP,the search ef ficiency of ALGA-SQP can be greatly enhanced.

(3)A hybrid multi-objective optimization algorithm,NSGA2-GAM,is presented,which solves the multi-objective optimization problem that takes into account both the interception time and interception uncertainty.Based on NSGA2,the algorithm expands the search space of NSGA2 by adding the optimal solutions of multiple single objective functions to the initial population,which makes the Pareto frontier searched by the NSGA2 has better boundary convergence and richer solution set diversity.At the same time,by using GAMto transform the multi-objective function into a mixedfunction for solving,the de ficiency of the local search ability of NSGA2 algorithm can be made up,and then a higher quality Pareto solution set can be obtained.

2.Problem description and analysis

Tangent impulse interception limits the impulse direction of the satellite.Fig.1 shows the schematic diagram of the entire tangent impulse interception process.

The satellite is initially located at P,and the target is initially located at P.After the satellite drifts from point Pto point P,a tangent impulse is applied to the interceptor at point P,so the intercepting orbit and the satellite orbit are tangent at point P.Let the transfer time of the satellite from the initial point Pto the impulse point Pbe t,the transfer time of the target from the initial point Pto the interception point Pbe t,and the transfer time of the interceptor from the initial point Pto the interception point Pbe t.If the interceptor wants to successfully intercept the target at P,it must satisfy the following equation:

Fig.1.Schematic diagram of the tangent impulse interception.

where?is the non-negative integer set and Tis the target orbital period.When equation(1)holds,ηis equal to the number of cycles of the target travels in its orbit.

From equation(1),we can see thatηis only a function of the satellite's true anomaly fat the impulse point Pand the target's true anomaly fat the interception point P.Therefore,for the minimum time interception problem,that is,solving a set of（f，f）that makes tminimum under the condition of satisfying equation(1).

The above-mentioned nominal orbit is designed to minimize the total interception time of the interceptor.Due to orbit determination error and impulse error,the actual flight path of the interceptor will deviate from the nominal orbit.So when we design the orbit,we also need to take the deviation factor into account.

After considering the orbit determination error,the actual state of the satellite at the impulse point can be expressed as:

After considering the impulse error,the actual velocity impulse obtained by the interceptor at the impulse point can be expressed as:

whereαandβare constant coef ficients.

When the satellite selects different impulse points,the error change of the interceptor from the impulse point to the interception point is shown in Fig.2.

It can be seen from Fig.2 that when the interceptor launches at different impulse points,terminal error distribution at the interception point will show a large difference.This is mainly caused by two factors.First,due to different impulse points correspond to different predetermined interception points,different interception orbits will lead to different transfer times for the interceptor to reach different interception points.Under the condition of the same initial error,a longer transfer time twill cause the interceptor's terminal error distribution to diverge more seriously at the predetermined interception point.Secondly,due to the different predetermined impulse speeds corresponding to different impulse points,a larger initial impulse speed will bring a greater initial speed error to the interceptor.Although a larger velocity impulse may shorten the overall transfer time,a larger initial velocity impulse will bring a larger initial velocity error,which may in turn make the terminal error distribution of the interceptor more divergent when it reaches the predetermined interception point.

From the above analysis,we know that the interceptor's terminal error distribution is a function of the interceptor's transfer time t,the impulse velocityΔv（t）,and the satellite's orbit determination errorδX at the impulse point.Therefore,to find the minimum terminal error of the interceptor is to solve a set of（f，f）,so thatσ（t）that meets the condition of Eq.(1)is minimized.If we want to take into account both the interception time and the interception uncertainty,we need to construct a multi-objective optimization function on（f，f）.By obtaining a series of optimal（f，f）,the transfer time tof the interceptor and its terminal error distributionσ（t）can reach Pareto Optimality.

Fig.2.Schematic diagram of the terminal position error distribution of the interceptor when launched at different impulse points.

3.Minimum time tangent impulse interception

Tangent impulse interception is the background of the problem in this paper,so in this section,we give a speci fic solution method for minimum time tangent impulse interception problem.We first use the flight path angle to express the interceptor's transfer time equation and give the existence range of the tangential impulse interception solution.When the impulse point is not fixed,we solve the global minimum time interception solution by optimizing the position of the initial impulse point and give an iterative calculation method of the minimum time tangent impulse solution.At the same time,to obtain the optimal solution more quickly and accurately,we also establish a global optimization model within the period of the satellite and give an optimization method based on the ALGA-SQP method to solve the minimum time interception solution.

3.1.Expression of interceptor transfer time equation based on the flight path angle

The traditional Lagrange transfer time equation has nothing to do with velocity direction.Therefore,to solve the tangent impulse interception problem,we first need to use the flight path angle to represent the transfer time equation.According to the literature

[26],by de fining the intermediate variableλ,we can express the interceptor's transfer time equation as the following three different forms according to the shape of the interceptor's orbit:

where Nrepresents the number of cycles the interceptor is running on the intercepting orbit.

where v（t）is the instantaneous velocity of the interceptor after the instantaneous impulse is applied;whenλ∈（0，2）,the intercepting orbit is an elliptical orbit;whenλ∈（2，∞）,the intercepting orbit is a hyperbolic orbit;whenλ=2,the intercepting orbit is a parabolic orbit.The expressions of the flight path angleγand the intermediate variableλare:

According to the change of the true anomaly of the satellite and target in the whole interception process,the transfer time of the satellite and target can be calculated:

where i=1 represents the satellite motion,i=2 represents the target motion;his the orbital momentum moment of the satellite or target;μis the gravitational constant;and fis the true anomaly of the satellite or target at the initial moment.Then,according to Eqs.(1),(7)and(8),we can calculate the tangent impulse interception solution at the given impulse point.

3.2.Tangent impulse interception solution with minimum time

Therefore,the interception time tof the interceptor is only a one-variable function of the impulse position f.Accordingly,we can solve the global minimum time interception solution by optimizing the initial impulse point position.The speci fic iterative calculation method is shown in Algorithm 1.

Although Algorithm 1 can be used to iteratively calculate the minimum time interception solution,the overall computation amount is too large.In order to get the global optimal impulse transmission scheme more easily and quickly,we establish an optimization model to solve the minimum time tangent impulse interception solution and use the optimization method to solve the problem.According to the previous description,we take the true anomaly fof the satellite at the impulse point and the true anomaly fof the target at the interception point as variables,take the impulse direction,impulse energy and the existence range of the interception solution as constraints,and take the minimum interception time tas the objective function,so as to establish the optimization model in a satellite cycle as follows:

Algorithm 1:Iterative calculation method of minimum time interception solution(1)Given the value range[f10，f1，max]of satellite impulse point f1 and iteration intervalΔf1.(2)For f1=f10:Δf1:f1，max do 1)According to Eqs.(9)and(10),calculate the existence interval（f2，s1，f2，s2）of the interception solution corresponding to the impulsive point f1.2)According to Eq.(11),calculate all the extreme points ofη,and divide the interval（f2，s1，f2，s2）according to the extreme points,so as to obtain all the solutions by the secant method.3)Delete all the solutions that do not satisfy the energy constraints and preserve the remaining solutions in the solution set.End(3)Compare the interception time t3 of all solutions in the solution set.The solution corresponding to the minimum t3 is the minimum time interception solution.

In order to find out the tangent impulse interception solution with the global minimum time accurately and ef ficiently,we use the hybrid optimization algorithm ALGA-SQP to solve the optimization model(12).The hybrid optimization algorithm combines the strong global search ability of the global optimization algorithm and the ef ficient local search ability of the local optimization algorithm.By taking the solution obtained by the ALGA algorithm as the initial value of the SQP search,the precise global optimal solution can be found quickly.AlGA deals with the nonlinear constraints by constructing sub-problems as shown in Eq.(13):

where f is the set of variables;ιis the estimated value of Lagrange multiplier;s is the displacement;κis the penalty coef ficient;Γ（f）is the objective function;cq（f）is the i-th inequality constraint;ceq（f）is the i-th equality constraint.

In this way,a new fitness function can be constructed by integrating the objective function and nonlinear constraints of the problem,and then GA can be used to solve the new fitness function.As for the SQP algorithm,there are many variations and extensions as a relatively mature nonlinear programming algorithm.In view of the fact that the optimization model established in this article is not unique,we only use the classical SQP algorithm[27]to solve the model and realize the algorithm solution by using the fmincon function that comes with Matlab.The complete algorithm process is shown in Algorithm 2:

Algorithm 2:The optimization method for finding the minimum time interception solution based on the ALGA-SQP algorithm(1)Given the initial population number,crossover and mutation probability,maximum genetic algebra,maximum impulse velocity increment,function tolerance,termination tolerance,and other initial parameters.(2)Randomly generate the impulse point time f1 in the range of[f10，f10+2π）and the corresponding interception time f2 in the range of（f2，s1，f2，s2）as the initial population.(3)While(genetic algebra≤maximum genetic algebra)&(average relative change of fitness function value≥function tolerance)1)According to the new fitness function(substituting Eq.(12)into Eq.(13)),calculate the fitness value of each individual in the population.2)Select,crossover and mutate the f1 and f2 codes of each group.3)Reorganize the code segments to get a new scheme,and update theι,s,andκaccording to the obtained f1 and f2,so as to modify the fitness function.End(4)Construct the QP sub-problem of Eq.(14).The optimal solution obtained by the ALGA is used as the initial search value to solve the QP sub-problem.Then we can get the search direction d k and the Lagrange multipliersαk andβk.(5)While(number of iterations≤maximum number of iterations)&(function value change≥termination tolerance)1)If d k=0 Output results,break End If 2)Use the golden section method to determine step size?k and let g k+1=g k+?k d k.3)Use the Davidon-Fletcher-Powell formula to modify the positive de finite approximation matrix H k of Hessen matrix of Lagrange function.End(6)Output the minimum time interception solution obtained by the SQP algorithm.

4.Optimal tangent impulse interception considering the terminal error

By using the method in the previous section,we can obtain the nominal orbit of minimum time interception under the condition of tangent impulse.In this section,we take the deviation factor into account in the whole orbit design process,so as to solve the orbit with the minimum error under the tangent impulse interception.Firstly,we establish the universal time equations for three different conic orbits.In this way,the universal variables at the interception point can be obtained through iterative calculation,and the interceptor state at the interception point can be obtained according to the Lagrange coef ficient.Then,in order to make the sampling result more accurate,we also give a method to calculate the terminal error distribution of the interceptor by using the bootstrap resampling technology.Finally,according to the geometric relationship between the error ellipsoid and the error box of the interceptor's terminal position,we designed an index to measure the distribution of the interceptor terminal error and established an optimization model to solve the minimum tangent impulse interception error.To get the global optimal solution quickly and accurately,the ALGA-SQP algorithm is still used to solve the optimization model.

4.1.Terminal error distribution calculation of the interceptor based on the universal time equation and bootstrap resampling

To calculate the terminal error of the interceptor at the interception point,it is necessary to adopt an appropriate deviation propagation method.Considering the in fluence of the nonlinearity,the analytical method based on STT is usually used to calculate the orbit deviation propagation.However,the trajectory shape of the interceptor may be an ellipse,hyperbola,parabola,and so on.In order to obtain the terminal error distribution of the interceptor more conveniently and accurately,we adopt the deviation propagation analysis method based on sample points.

The dynamic models of the three kinds of conic orbits are different.In order to unify the trajectory extrapolation models of the interceptors,we use the universal variableχinstead of the time variable to unify the motion equations of interceptors with different orbital forms.

To avoid the singularity caused by the semi-major axis of the parabola orbit being in finite,we de fine the reciprocal of the semimajor axis of the interceptor orbit as follows:

where,if a＞0,the interceptor trajectory is elliptical;if a=0,the interceptor trajectory is parabolic;if a＜0,the interceptor trajectory is hyperbolic.

According to the Ref.[28],the universal time equation of the interceptor orbit can be described as:

where ois the intermediate variable,o=r·v/μ;rand vare the position and speed of the interceptor at the launching time;ris the interceptor position at the end of orbit extrapolation;U,U,Uand Uare universal functions,and their expressions are:

The Lagrange coef ficients expressed by the universal function are:

where G,G,Gand Gare Lagrange coef ficients.

Then substituting the universal variableχobtained by the iterative calculation into Eq.(17),the position and velocity of the interceptor at the extrapolation end point are as follows:

To sum up,we calculate the terminal error distribution of the interceptor according to the following three steps:(1)generate a set of initial deviation samples of the interceptor according to Eqs.(2)-(4);(2)obtain the state of the sample point at the extrapolation endpoint through the universal time equation;(3)count the first two order statistical moment information of all the terminal state sample points,that is,the mean value︿mand the variance^Cof the position state error of the terminal state sample points.

In order to ensure the accuracy of the final statistical information,a large number of sample points are needed in the calculation process.However,if the sample size is too large,it will consume a lot of computing resources.Therefore we use the bootstrap resampling method to reduce the estimation error of each statistical moment.In this way,the accuracy of the low-order statistical information can be guaranteed as much as possible under the premise of using fewer samples.The bootstrap method is a resampling method with putting back[29].This method does not make any assumptions about the overall distribution of the unknown samples,but does resampling of the existing samples to generate a large number of regenerated samples,so as to simulate the unknown distribution.At the same time,the same operator is used to obtain the estimator of the regenerated sample,and then further statistical work is carried out based on the estimator.The statistical information calculation method of the interceptor terminal position state error based on the bootstrap:(1)obtain the state of all sample points at the extrapolation terminal point by the universal time equation;(2)form all the state sample points at the extrapolation endpoint into a new set of samples,and obtain B groups of bootstrap samples by sampling with the return;(3)calculate the first two statistical moments︿mand^Cof the position error of each Bootstrap sample group;(4)take the mean value of all statistical moments of the B groups as the final position state error distribution estimation of the interceptor:

4.2.Optimal interception solution of tangent impulse considering the terminal error

To optimize the interception accuracy of the interceptor at the interception point,we need to design an index to evaluate the terminal error of the interceptor.At the same time,for the interception problem,compared with the terminal velocity error of the interceptor,we should pay more attention to whether the interceptor can accurately hit the target at the interception point.Therefore,we choose the standard deviation of the interceptor's position error at the interception point to measure the terminal error distribution of the interceptor.

In the Ref.[30],Gang Zhang and Yunhai Geng used the terminal stand deviation magnitude to indicate the terminal stand deviation magnitude.To make the quantitative indicators of uncertainty more physical,we chose the volume of error ellipsoid to describe the terminal deviation visually.Taking the projection of the error on a two-dimensional plane as an example,we show the geometric relationship between the standard deviation of the interceptor position error and the error ellipse in Fig.3.It can be seen from Fig.3 that 2σ（l=x，y）geometrically represents the edge length of the circumscribed rectangle of the error ellipse.Similarly,2σ（l=x，y，z）represents the side length of the circumscribed cuboid of the error ellipsoid geometrically.

Fig.3.Projection of the error on a two-dimensional plane.

When optimizing the terminal error distribution of the interceptor,we can indirectly minimize the volume of the error ellipsoid by minimizing the volume of the circumscribed cuboid.In this way,it is not necessary to calculate the volume of the error ellipsoid directly,thus avoiding the eigenvalue decomposition of the error covariance matrix.Therefore,we design the optimization index of the interceptor terminal position error distribution as the volume of the circumscribed cuboid of the position error ellipsoid:

In this way,we take the true anomaly fof the satellite at the impulse point and the true anomaly fof the target at the interception point as variables,take the impulse direction,impulse energy and the existence range of the interception solution as the constraints,and take the minimum volume Vof the circumscribed cuboid of the interceptor position error ellipsoid as the objective function,so as to establish the optimization model for solving the minimum tangent impulse interception error in a satellite cycle as follows:

whereσ（l=x，y，z）is calculated by Eq.(19).

To make the optimization process ef ficient and accurate,we still use the hybrid optimization algorithm ALGA-SQP to solve the optimization model.By combining the objective function and nonlinear constraints in equation (21),we can solve the constraint problem in the optimization process.Then the GA algorithm with a strong global search ability is used to find the relatively optimal global solution of the new fitness function.By taking the solution as the initial value of the SQP search,a more accurate global optimal solution can be quickly found.The complete algorithm flow is shown in Algorithm 3.In order to make the description more rigorous,it is hereby explained:the solution obtained by the ALGA-SQP algorithm is not strictly a globally optimal solution,but a numerical solution that's very close to the optimal global solution.

Algorithm 3:The optimization method for minimizing the interceptor terminal position error based on the ALGA-SQP algorithm(1)Given the initial population number,crossover and mutation probability,maximum genetic algebra,maximum impulse speed increment,the standard deviation of the satellite orbit determination error at the impulse pointσδX（tP1）,coef ficient of the standard deviation of the impulse velocity errorα，β,Monte Carlo shooting times,bootstrap resampling times B,function tolerance,terminal tolerance,and other initial parameters.(2)Randomly generate the impulse point time f1 in the range of[f10，f10+2π）and the corresponding interception time f2 in the range of（f2，s1，f2，s2）as the initial population.(3)While(genetic algebra≤maximum genetic algebra)&(average relative change of fitness function value≥function tolerance)1)For each individual in the randomly generated initial population,according to the satellite orbit determination error standard deviationσδX（tP1）and the impulse velocity error standard deviation coef ficientα，βat the impulse point,randomly generate a group of interceptor deviation samples at the impulse point by equations(2)-(4).2)According to the universal time equation,calculate the terminal position error distributions of the interceptor in three axial directions by the Monte Carlo shooting method.3)Use the Bootstrap resampling method to calculate the final standard deviation estimator of the interceptor state error and the final volume estimation Vbox of the circumscribed cuboid of the interceptor position error ellipsoid.4)According to the new fitness function(substituting equation(21)into equation(13)),calculate the fitness value of each individual in the population.5)Select,crossover,and mutate the f1 and f2 codes of each group.6)Reorganize the code segments to get a new scheme,and update theι,s,andκaccording to the obtained f1 and f2,so as to modify the fitness function.End(4)Take the optimal solution obtained by the ALGA as the initial search value,and use the SQP to solve equation(21)(the speci fic solving method is consistent with Algorithm 2).(5)Output the minimum terminal position error solution of the interceptor obtained by the SQP precise search.

Fig.4.Geometric representation of the goal attainment method.

5.Pareto optimal solution set considering interception time and terminal uncertainty under tangent impulse interception

When solving the minimum interception time t,the optimal solution obtained by the optimization model in Eq.(12)tends to select the impulse point fwhich needs a higher impulse velocity,so as to reduce the overall interception time of the interceptor.From Eq.(4),it can be seen that the impulse velocity error of the satellite is directly proportional to the impulse velocity of the satellite.The optimal impulse point obtained by optimizing the terminal error of the interceptor according to formula(21)may not make the interceptor have the shortest interception time.Therefore,considering the interception time and the interception uncertainty,we take the interception time tand the volume Vof the circumscribed cuboid volume of the interceptor terminal position error ellipsoid as the objective function and establish the following multi-objective optimization model:

By Eq.(22),we transform the tangent impulse interception problem into a multi-objective optimization problem.Different from single-objective optimization,the improvement of some objectives in multi-objective optimization may cause the deterioration of others.In this regard,the optimal solution is usually expressed by the non-inferior solution.Multi-objective optimization requires the algorithm to find as many and evenly distributed solutions as possible in the set of non-inferior solutions,that is,the Pareto front.Compared with other mathematical programming methods,the genetic algorithm(GA)can get multiple solutions in parallel after one operation because of its group operation,so it is very suitable for solving multi-objective optimization problems.Moreover,GA is not sensitive to the shape and continuity of the Pareto front-end and can approach the optimal front-end without convexity or discontinuity.Therefore,we use the improved NSGA2-GAM hybrid optimization algorithm to solve Eq.(22).

NSGA2 algorithm uses the non-dominated solutions to sort quickly,and uses the elite strategy and crowding degree screening to obtain the Pareto optimal solution set.At the same time,the NSGA2 algorithm does not need to provide the initial value,which can avoid the subjectivity of the coef ficient weighting method.However,the simulated binary crossover operator used in the NSGA2 algorithm is limited in search range and prone to local optimization.To solve this problem,we use the solution of the single-objective function to start the NSGA2 algorithm.In other words,when NSGA2 is initialized,the initial population containing the optimal solutions of Eq.(12)and Eq.(21)is generated,so as to expand the search space.Meanwhile,in order to further enhance the accuracy of the solution set,after obtaining the Pareto front by the NSGA2,we use the goal attainment method to transform the multi-objective optimization function into a mixed function.Then,the SQP algorithm is used to solve the mixed-function,and the Pareto optimal solution set with higher quality is obtained.In this way,the hybrid optimization algorithm can quickly search for Pareto frontier with high quality by combining the global optimization algorithm NSGA2 with strong global searchability and the local optimization algorithm SQP with high local searchability.

The goal attainment method transforms the multi-objective optimization function into:

In addition,when using the NSGA2 algorithm to deal with the constraints in Eq.(22),we use the Penalty algorithm to replace the augmented Lagrange method.If the individual satis fies the constraints,the penalty function is the fitness function.If the individual does not satisfy the constraints,the penalty function is the sum of the maximum fitness functions of the feasible members and the constraint violation of the infeasible individuals.

To sum up,Algorithm 4 gives the complete process of solving the tangent impulse interception problem by the NSGA2-GAM algorithm.

Algorithm 4:NSGA2-GAM multi-objective optimization method for solving the tangent impulse interception problem(1)Given the initial population number,crossover and mutation probability,maximum genetic algebra,maximum impulse speed increment,the standard deviation of the satellite orbit determination error at the impulse pointσδX（tP1）,coef ficient of the standard deviation of the impulse velocity errorα，β,Monte Carlo shooting times,bootstrap resampling times B,maximum stall generations,function tolerance,terminal tolerance,and other initial parameters.(2)Randomly generate the impulse point time f1 in the range of[f10，f10+2π）and the corresponding interception time f2 in the range of（f2，s1，f2，s2）.Take the generated series(f1,f2)and the minimum time solution and minimum terminal error solution obtained by optimization as the initial population.(3)While(genetic algebra≤maximum genetic algebra)&&(geometric mean value of the relative variation of the Pareto solution spread after exceeding the maximum stall generations≥function tolerance｜｜final spread≥average spread after exceeding the maximum stall generations)1)Use Penalty algorithm to calculate the virtual fitness of the individuals,and use the fast non-dominated sorting method to rapidly stratify the population;2)Based on the two attributes of each individual's non-dominance level and crowdedness,all individuals in the population are screened by the tournament selection method to obtain the parent population 3)The elitists in the parent population and the offspring population are mixed to form the next generation population.End(4)Construct a mixed-function according to the Pareto frontier obtained by the NSGA2,and calculate the target value Fk*（f p）and corresponding weight w*k*（f p）for each non-inferior solution.Then transform the interception problem into a single-objective optimization problem through the goal attainment method,and solve the problem by SQP.(5)Output the optimal Pareto solution set of tangent impulse interception obtained by the hybrid optimization algorithm.

6.Analysis of examples

To verify the effectiveness of the optimization method proposed in this paper,we carry out the simulation veri fication in three parts.They are simulation studies on the minimum time tangent impulse interception without considering the interception error,the optimal tangent impulse interception only considering the terminal error,and the multi-objective optimization of the tangent impulse interception considering both the interception time and the terminal uncertainty.

In the example,the upper limit of the satellite impulse velocity is 3 km/s.The initial orbit elements of the satellite and target are given in Table 1.

Table 1 Initial orbital elements of the satellite and target.

6.1.Simulation studies on the minimum time tangent impulse interception without considering the interception error

According to Algorithm 2,we use the ALGA-SQP hybrid optimization algorithm to search for the minimum time interception solution which satis fies the constraint requirements of equation(12)in a satellite cycle.Set the initial population number as 30,crossover probability as 0.80,mutation probability as 0.20,maximum genetic algebra as 30,maximum impulse velocity increment as 3 km/s,initial penalty as 10,penalty factor as 100,function tolerance as 10and termination tolerance as 10.The iteration process is shown in Fig.5.

To test the accuracy of the global hybrid optimization algorithm,we iteratively calculate all the feasible solutions in a satellite cycle according to Algorithm 1.Take the iteration interval as 1.When there are multiple solutions at the same impulse point,we take the minimum time interception solution satisfying the constraint requirements.In this way,each fixed impulse point f∈[f，f+2π）in the satellite cycle corresponds to only one corresponding terminal interception point f.The iterative calculation results are shown in Fig.6.

Fig.5.Iterative processes of solving the minimum time tangent impulse interception solution by the ALGA-SQP algorithm.

As can be seen from Fig.6,there are two breakpoints for the interception solution curve within a satellite cycle.The first breakpoint is located at f=429.From here on,the number of running cycles Ncorresponding to the minimum time interception solutions satisfying the constraint requirements is changed from 1 to 0.The second breakpoint is at f=463.From this point,the trajectory of the interceptor corresponding to the minimum time interception solutions satisfying the constraint requirements change from hyperbola to ellipse.From Fig.6(a),we can see that the satellite initial impulse point corresponding to the minimum time interception solution satisfying the constraint requirements is f∈（462，463）.Therefore,it can be seen that the hybrid optimization solution obtained by the ALGA-SQP is very close to the global optimal solution,which veri fies the effectiveness of the method in this paper to solve the minimum time tangent impulse interception solution.

Besides,to verify the superiority of the ALGA-SQP algorithm,we will compare computing time for the iterative calculation method and the ALGA-SQP algorithm.In order to save calculation time,we choose the method in the appendix of the literature[7]when calculating the transfer time of satellite and target.In the hardware environment of Intel Core i7-8750H CPU and 16 GB memory,using the Matlab 2019b software platform,the total calculation time consumed by the above iterative calculation process is 753.16s.In the same hardware and software environment,the total calculation time of the ALGA-SQP algorithm is only 1.92s.It is mainly because all feasible solutions in a satellite period in sequence during the iterative calculation need to be calculated.Since the selected time interval is 1,the complete iterative operation needs to be repeated 361 times.Furthermore,we need to use the secant method to solve in each iteration,so the total calculation time of the iterative method is longer.In contrast,ALGA-SQP has completed convergence after only seven iterations in total.And parallel calculations are used throughout the calculation process,so the time spent is shorter.Compared with iterative calculation,using the ALGA-SQP hybrid optimization algorithm can signi ficantly reduce the calculation time while ensuring the calculation accuracy.

The whole interception process corresponding to the minimum interception time solution obtained by the hybrid optimization is shown in Fig.7.

Fig.6.Iterative solutions within a satellite cycle.

Fig.7.Minimum time interception process.

6.2.Simulation study on the optimal tangent impulse interception problem considering terminal error only

When the satellite chooses different impulse points,we can get the state error distribution of the interceptor at the terminal position according to the universal time equation.The calculated interceptor state errors corresponding to different impulse points in a satellite cycle are shown in Table 2:

It can be seen from Table 2 that when the sampling number is 1000,comparing with the calculation results with the sampling number of 20,000,the bias of the position standard deviation of the interceptor in the three-axis directions is 1.05%,2.03%,and 1.74%respectively,and the bias of the circumscribed cuboid volume of the ellipsoid is 2.71%.From this,we can see that the maximum bias of the standard deviation of the terminal position error of the interceptor and the maximum bias of the volume of the circumscribed cuboid of the interceptor position error ellipsoid are less than 3%under different impulse points in a satellite period.So in the next simulation,we choose the sample number as 1000.

Table 2 Terminal position error of the interceptor under different impulse points.

Table 3 Optimal Pareto solution set obtained by the GAM optimization algorithm.

According to Algorithm 3,we use the ALGA-SQP hybrid optimization algorithm to search for the minimum interceptor terminal position error solution which satis fies the constraint requirements of Eq.(21)in a satellite cycle.Set the initial population number as 30,crossover probability as 0.80,mutation probability as 0.20,maximum genetic algebra as 30,maximum impulse velocity increment as 3 km/s,Monte Carlo shooting times as 1000,resampling number B of the Bootstrap as 50,initial penalty as 10,penalty factor as 100,function tolerance as 10and termination tolerance as 10.The iteration process is shown in Fig.8.

It can be seen from Fig.8 that ALGA converges after five iterations,and the minimum interceptor terminal position error solution obtained by searching in the global range is（f，f）=（457.8950，168.8996）.The standard deviation of the terminal position error of the interceptor in three coordinate axes isσ=6.6275km,σ=0.5291km andσ=2.6639km,respectively;the volume of the circumscribed cuboid of the interceptor position error ellipsoid is V=74.7240 km.The transfer time tof the interceptor is 6679.83s,and the impulse energy isΔv（t）=1.4494km/s;the total transfer time tof the target is 13868.56s,λis 1.8279,and the intercepting orbit shape is ellipse;the corresponding maximum constraint violation is 2.0522×10.Then we take this solution as the initial search value of the SQP algorithm.After four iterations,a hybrid optimization solution with smaller terminal position error is obtained near the solution,and the optimal solution is（f，f）=（458.0244，168.2538）.The corresponding standard deviation of the terminal position error of the interceptor in three coordinate axes isσ=6.3826km,σ=0.5129km andσ=2.6410km,respectively;the volume of the circumscribed cuboid of the interceptor position error ellipsoid is V=69.1663 km.The transfer time tof the interceptor is 6520.58s,and the impulse energy isΔv（t）=1.5km/s;the total transfer time tof the target is 13685.64s,λis 1.8516,and the intercepting orbit shape is ellipse;the corresponding maximum constraint violation is 3.8765×10.

Fig.8.The iterative process of solving the minimum interceptor terminal position error solution by the ALGA-SQP algorithm.

To verify the accuracy of the global hybrid optimization algorithm,we iteratively calculate the feasible solutions with the minimum terminal position error corresponding to all impulse points in a satellite cycle.Take the iteration interval as 1,Monte Carlo shooting times as 1000,resampling number B of the Bootstrap as 50.When there are multiple solutions corresponding to the same impulse point,we take the solution which has the minimum terminal interception error and satis fies the constraint requirements. In this way, each fixed impulse point f∈[f，f+2π）in the satellite cycle corresponds to only one corresponding terminal interception error.The iterative calculation results are shown in Fig.9.

As can be seen from Fig.9,there are three break points in the minimum terminal interception error solution curve within a satellite cycle.The first breakpoint is located at f=429.From here on,the number of running cycles Ncorresponding to the minimum terminal interception error solutions satisfying the constraint requirements is changed from 1 to 0.The second breakpoint is located at f=458.From this point,the standard deviation coefficientαof impulse velocity error corresponding to the minimum terminal interception error solution satisfying the constraint requirements is changed from 0.0005 to 0.001.The third breakpoint is located at f=463.From here on,the orbit of the interceptor corresponding to the minimum terminal interception error solution satisfying the constraint changes from the hyperbola to the ellipse.From Fig.9(d),we can see that the satellite initial impulse point corresponding to the minimum terminal interception error solution satisfying the constraint requirements is f∈（458，459）.Therefore,it can be seen that the hybrid optimization solution obtained by the ALGA-SQP is very close to the global optimal solution,which veri fies the effectiveness of the method in this paper to solve the minimum terminal interception error solution.

The minimum terminal position error distribution of the interceptor after hybrid optimization is shown in Fig.10.

Fig.9.Iterative calculation solutions of the interceptor terminal position error in one satellite cycle.

6.3.Simulation study on the multi-objective optimization of the tangent impulse interception considering both the interception time and the terminal uncertainty

By comparing the simulation results of the ALGA-SQP hybrid optimization algorithm,we can see that the interception time corresponding to the global minimum time interception solution（462.6463，160.1618）is 2410.13s faster than that of the global minimum terminal position error solution（458.0244，168.2538）.And the larger impulse velocity error also makes the terminal error distribution more divergent when the interceptor reaches the predetermined interception point.

Therefore,to balance the interception time and the terminal position error,we use the NSGA2-GAM hybrid optimization algorithm according to algorithm 4 to search a Pareto optimal solution set which satis fies the constraint requirements of Eq.(22)in a satellite cycle.Set the initial population number as 150,crossover probability as 0.80,mutation probability as 0.20,maximum genetic algebra as 200,maximum impulse velocity increment as 3 km/s,Monte Carlo shooting times as 1000,resampling number B of the Bootstrap as 50,constraint tolerance and function tolerance of NSGA2 and GAMboth as 10and 10.The optimization results of the NSGA2 are shown in Fig.11.

As can be seen from Fig.11,the NSGA2 algorithm converges after 122 iterations.The average spread is 0.7638 at the end of the algorithm.We can see that the Pareto solution set is mainly concentrated near the minimum time solution（462.6463，160.1618）and the minimum terminal position error solution（458.0244，168.2538）.It can be seen that for the simulation examples in this article,the impact of the interceptor transfer time on the interceptor terminal position error is far greater than the impact of the impulse velocity error.Therefore,the Pareto solution set obtained in the example shows a discrete step form.Near the minimum time solution,due to the short interception time,the impulse velocity error plays a leading role in affecting the interceptor terminal position error,so some Pareto solutions appear here.In the vicinity of the minimum terminal position error solution,the impulse velocity reaches 1.5 km/s,which causes the standard deviation coef ficient of the impulse velocity error to increase,and the in fluence of the impulse velocity error to increase,so some Pareto solutions also appear here.Therefore,the Pareto solution set obtained by simulation is only concentrated in the vicinity of the minimum time solution and the minimum terminal position error solution.

To get better optimization results,we add the minimum time solutions and the minimum terminal position error solutions to the initial population of the NSGA2 optimization algorithm,so as to expand the search space and further enhance the accuracy of the solution set.The optimization results after modifying the initial population are shown in Fig.12.

Fig.10.The minimum terminal position error distribution of the interceptor.

Fig.11.Optimization results of the NSGA2.

Fig.12.Optimization results of the NSGA2 after modifying the initial population.

Fig.13.Optimization results of the GAM.

Fig.14.Optimal compromise solutions.

It can be seen from Fig.12,the NSGA2 algorithm converges after 169 iterations.The average spread is 0.9788 at the end of the algorithm.It can be seen from Fig.12 that when the optimal solution of the single objective function is added to the initial population,compared with the optimization results of the randomly initialized population,the global Pareto frontier obtained by NSGA2 optimization algorithm has better boundary convergence and richer solution set diversity.Moreover,the Pareto solutions in the solution set have shorter overall interception time and smaller terminal position error.The results show that if the algorithm starts from the optimal solution of the single objective function,it can make up for the defects of NSGA2,such as low search range and easy to be limited to local optimization.

Fig.15.Comparison between the interception time compromise solution and the interception error compromise solution.

To complement the NSGA2 algorithm,further enhance the local search ability of the algorithm.Next,we use the GAMalgorithm to transform the multi-objective function into a mixed-function,take the global search results as the initial values,and use the SQP algorithm to further search for the solutions of the mixed-function,so as to obtain the Pareto solution set with higher quality.The weight of each non-inferior solution can be calculated as{1.0320,1.0214,0.9996,0.9861,1.0039,1.0111,0.9662,0.9788,0.9680,0.9697}.The weights of the two objective functions at each noninferior solution are {(0.0337,1.0313),(0.0224,1.0210),(0,0.9992),(0.9723,0),(0.0043,1.0034), (0.0117,1.0107),(0.9336,0),(0.9581,0),(0.9369,0),(0.9404,0)}。Therefore,at each non-inferior solution,the multi-objective optimization function can be combined into a mixed optimization function according to the above weight values.Take the optimization results as the initial search values and the mixed-function at each non-inferior solution as the objective function.With the constraints unchanged,the SQP algorithm is used to further optimize each non-inferior solution.The hybrid optimization results are shown in Fig.13.

The obtained optimal Pareto solution set is shown in Table 3.

It can be seen from Fig.13 and Table 3 that the Pareto frontier obtained by the GAM has better boundary convergence compared with the optimization results obtained by the NSGA2.Meanwhile,due to the increase of the function tolerance,compared with the global optimal solution obtained by the single-objective optimization function,there are solutions with shorter total interception time and smaller terminal position error in the Pareto solution set obtained by the NSGA2-GAM optimization algorithm.The results show that the hybrid optimization algorithm NSGA2-GAM can overcome the defect of weak local search ability of the NSGA2 and obtain a better Pareto frontier.

By comparing all the solutions in the Pareto solution set,it can be seen that although a larger impulse velocity makes the interceptor's overall transfer time shorter.However,the more signi ficant initial velocity error caused will make the interceptor's terminal error distribution more divergent.Therefore,when interception requires a smaller terminal position error,we choose the Pareto solution with the smallest terminal position error as the interception error compromise solution.This solution is also the first endpoint of the Pareto front.Its overall interception time tis 6518.43685s,and its terminal position error Vis 66.80674 km.Comparing the solutions in the lower right part of Fig.13,we can see that increasing the terminal position error will basically not reduce the interceptor's overall interception time.Compared with the minimum time Pareto solution,the leftmost solution in the lower right part of Fig.13 has a smaller terminal interception error under the premise of basically the same interception time.Its overall interception time tis 4107.89010s,and the terminal position error Vis 136.29004 km.Compared with the minimum time Pareto solution(another endpoint of the Pareto front,the rightmost solution in the lower right part of Fig.13),this solution's terminal position error is shortened by 0.35 km,and the overall interception time is only increased by 4.27×10s.The solution's terminal error is shortened by 0.26%,and the interception time is increased by 0.0104‰.For most calculation examples,the Pareto solution set's endpoints are usually selected as the optimal compromise solutions.However,in this example,the pros and cons of the two solutions are not very obvious.Therefore,based on actual magnitude,we choose the leftmost solution in the lower right part of Fig.13 as the interception time compromise solution.Compared with the interception error compromise solution,the terminal position error of the interception time compromise solution is increased by 69.48 km,but the overall interception time is shortened by 2410.55s.

The two optimal compromise solutions and their comparison results are shown in Fig.14 and Fig.15,respectively:

7.Conclusions

By taking the tangent impulse interception mission as the background,this paper starts by solving single-objective optimization problems such as the minimum time interception orbit and the minimum terminal interception error orbit and gives an orbit optimization method considering the interception time and the interception uncertainty.

The simulation example shows that:

(1)Compared with the iterative method,the ALGA-SQP hybrid optimization algorithm combines the advantages of ALGA and SQP.It can quickly search out the optimal global solution without providing the initial value.Therefore,it is more suitable for solving the minimum time tangent impulse interception problem and the minimum terminal interception error problem.

(2)The universal time equation uni fies three different conic orbits by using the universal variables,which effectively reduces the dif ficulty of solving the interceptor's terminal error.At the same time,Bootstrap resampling can make the statistical results of the terminal error more stable,thus reducing the number of Monte Carlo sampling and saving the cost of calculation.

(3)NSGA2 can search out the Pareto frontier satisfying the restrictions without providing the initial value.By adding the optimal solution of the single-objective function into the initial population,NSGA2 can search out the Pareto frontier with better boundary convergence and richer solution set diversity.Taking the search results of the NSGA2 as the initial value of the GAMalgorithm,the multi-objective function can be transformed into a mixed-function to solve,which further enhances the local search ability of the NSGA2 algorithm.Thus the Pareto solution set with higher quality is obtained.

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