An Adaptive Remeshing Procedure for Proximity Maneuvering Spacecraft Formations

We consider a methodology to optimally obtain reconfigurations of spacecraft formations. It is based on the discretization of the time interval in subintervals called the mesh and the obtainment of local solutions on them as a result of a variational method. Applied to a libration point orbit scenario, in this work we focus on how to find optimal meshes using an adaptive remeshing procedure and on the determination of the parameter that governs it.


Introduction
Formation flying concepts have been growing in interest during the last years.In different scenarios, formations or clusters of small satellites can perform like a virtual larger telescope, obtaining equivalent information, but with a reduced cost.Mission concepts like the NASA Terrestrial Planet Finder 1 , the ESA Darwin 2 , the NASA MAXIM 3 , or the ESA XEUS 4 are just few examples that remark the importance of this new technology for space telescopes.
Nevertheless, formation flying still demands many new technologies to be successfully implemented.Usually the spacecraft must be located and maintained within a very narrow range of relative distances, and severe constraints like this, increasing the already high complexity of the mission design see for instance the works of Farrar et al. 5 and references therein .Also there are many other issues that need to be addressed as well.For instance, a main one is collision avoidance when maneuvering or reconfiguring the formation.Typically, from one observation to the next one, the formation needs to change the pointing goal and eventually change its pattern.To this end, some representative techniques considered are

The FEFF Methodology
In this section we present a brief summary of the basics of the FEFF methodology that can be found fully developed in 9, 10 .This methodology was made with the purpose of systematically computing reconfigurations of spacecraft formations located in libration point orbits.However, it could also be generalized for formations about the Earth or for in free space.
As it is well known, the vicinity of the Lagrangian points L 1 and L 2 is a very convenient place for space observatories L 1 for the Sun and L 2 for deep space .In this paper we consider a formation of spacecraft located in a halo orbit about L 2 .We assume that the formation is made of N spacecraft flying with a particular pattern and our objective is to perform a reconfiguration in a fixed time T .We also assume that the spacecraft are in a small formation.This is, the distance between them is only of the order of few hundreds of meters, both in the initial and the final configurations.The objective of the FEFF methodology is to find a suitable trajectory for each of the spacecraft that delivers it to the goal position, with minimum fuel consumption and avoiding collisions with other spacecraft.
Since the formations are small with respect to the amplitude of the halo orbit, we consider the linearized equations of motion about the nonlinear orbit.In 11 we have studied the impact of the nonlinear part, concluding that, for orbits with a diameter of a few hundreds of meters i.e., the usual length for a formation , these linearized equations model the dynamics of the formation in a very good way.
Then, according to these hypotheses, associated to each spacecraft in the formation, we have an equation of the form where A t is a 6 × 6 matrix and X refers to the state of the spacecraft see The appendix .
Usually the origin of the reference frame for the X coordinates is the nominal point on the base halo orbit at time t being the orientation of the coordinate axis parallel to the ones of the RTBP model.
In order to perform the reconfigurations, we consider a control function for each of the spacecraft and also we include the boundary conditions corresponding to the initial and final states:

2.2
Here X 0 i and X T i stand for the initial and final states of the ith spacecraft inside the formation, and U 1 , . . ., U N , are the controls we are searching with the aim of being optimal in terms of fuel consumption.
The basis of the FEFF methodology is to use the finite element method in time to discretize the spacecraft trajectories and to obtain the controls.Essentially, for each spacecraft, the time interval 0, T which we consider for the reconfiguration is split in M i subintervals of the domain called elements.For a given trajectory, its mesh can have elements of different lengths and of course different satellites can have associated different meshes.This will depend on the nature of the trajectories of the reconfiguration and the path they follow.Using this mesh we impose that controls are maneuvers applied at the points where two elements join the nodes .The finite element methodology is used to formulate the problem and to obtain, by means of this discretization, a relation between the states of each spacecraft at the nodes of the elements and the Δv maneuvers applied.We note that in our discretization we use elements with two nodes located at the ends of each element consecutive elements share the connecting node .These elements are called linear elements and their associated truncation error is according to the linear approximation taken about the nonlinear orbit when considering 2.1 .Usually, in the finite element methodology, the solution inside each element is approximated by a polynomial and, for linear elements, this polynomial is of degree one.Finite elements in time and greater orders have also been considered in more general formulations of optimal control problems.An interesting presentation can be found in 12 .
By means of the procedure FEFF, we reduce the reconfiguration problem to an optimal control problem with constrains.The functional to be minimized is related to the fuel consumption and is taken as the sum of the norm of the maneuvers: where || * || denotes the Euclidean norm, v i,k is the maneuver applied at the kth node of trajectory, and i and ρ i,k are weighting parameters that can be used, for instance, to penalize fuel consumptions of selected spacecraft with the purpose of balancing fuel resources.For clarity, in 2.3 we consider that ρ i,k multiplies the modulus of the maneuver, but in a similar way we can impose a weight on each of its components.An important issue in the formulation of the procedure is collision avoidance between satellites.It enters in the optimization problem as constraints.We consider that each spacecraft is surrounded by a security sphere that cannot intersect during the reconfiguration process.Again, the discretization of the time interval made by the finite element methodology provides an efficient implementation to check these constraints.

Adaptive Remeshing Applied to Reconfigurations
The objective of this paper is how to systematically obtain good meshes for the reconfiguration problem.We note that, for a given mesh, the FEFF methodology computes the trajectory of the spacecraft in such a way that J 1 in 2.3 is minimized.But at the end, the trajectories of the satellites have been obtained after a discretization process and the error of the discretized approximated trajectories with respect to the exact solutions of the problem is not obvious.If we take a small number of elements, we can have a poor model that is not accurate enough.On the other side, as it is well known in the finite element methodology, the approximated trajectories converge towards the true solution when the diameter of the mesh the length of the longest time interval tends to zero.
Of course when we increase the number of elements in the mesh, it also increases the complexity of the computations, the required CPU time, and the representation of the solution someway.Also we could end up with ill conditioned problems when the number of elements is very high due to the presence of very small maneuvers.To overcome these difficulties, adaptive remeshing is a technique that allows us to work in the other way round.Fixing an acceptable level of accuracy, and by means of an iterative procedure, one can find "the coarsest mesh" providing approximate trajectories with the required accuracy.
Another issue we have to deal is related to the minimization of the functional 2.3 .Its derivatives are ill conditioned when one or more delta-v are near zero.Since the objective is to find these maneuvers as small as possible, one may expect problems if we perform the computations in a naive way.
We address these two facts in a two-step methodology.First we find an initial approximation of the solution minimizing an alternative functional.In a natural way we have chosen the functional: because it is also directly related to fuel consumption and it does not have any ill-conditioning problems when computing derivatives near zero.Using this target functional, there are no problems in finding a solution; moreover, the errors associated to a coarse mesh are not critical for the second step.Let us call FEFF-DV2 the procedure that provides this approximated solution.Starting with FEFF-DV2, we usually consider a mesh with a small number of elements generally from 5 to 10 and we take all of them of the same length.We note that since FEFF-DV2 is an optimization problem, we need an initial seed.For this purpose we consider each spacecraft alone, this is, we solve N independent problems, where we compute the trajectory which minimizes J 2 without taking into account possible collision risks.
Using the same discretization as in FEFF formulation, these initial seeds can be found semianalytically by means of solving a linear system 10, The proof .The solution is unique considering that the elements are all of the same length.Moreover, if the obtained trajectories do not have collision risks the exclusion spheres do not intersect , they are already the output of FEFF-DV2 i.e., the approximate solution for the given mesh that minimizes J 2 .
In the second step of the procedure, we consider an adaptive remeshing strategy with two purposes: to control the error due to the finite element methodology and to suppress all the nodes that can give ill-conditioning problems in the minimization of J 1 .Let us call FEFF-DV1 the formulation that uses the FEFF methodology to optimize the functional J 1 once the nodes that could give ill-conditioned problems have been removed.Then the general idea of the adaptive methodology follows the scheme displayed in Figure 1.Once we have the approximation given by FEFF-DV2, we start the iterative second step which involves FEFF-DV1 and an estimation of the error that is used to generate a new mesh.This iterative procedure is repeated till the estimated error is below the given threshold requirement.
Finally let us comment more in detail how the adaptive remeshing works.The general idea is that, given a threshold value e, we want to find a mesh that provides an approximate solution with error understood as the difference between the solution of the problem and its approximation inside of an element less than e in some norm.
Adaptive remeshing methods penalize the elements where the error is considered big, dividing them into smaller elements.On the other hand, if the estimation of the error is small in an element, this element is made bigger in the next iteration.Since, as we will see, our estimation of the error is basically related to the value of the delta-v maneuvers to be implemented; roughly speaking, this method tends to increase the length of the elements which have associated small delta-v and tends to decrease the length of the elements which have associated big delta-v's.As a consequence, it is also suitable to avoid the ill-conditioned problems that FEFF-DV1 might have.
Essentially, to decide whether the current mesh is good enough or not, we base on a criterion which compares the modulus of the estimated error, ||e||, on the mesh with the total gradient of the solution related in our problem with velocities .For this purpose we compute the following integral by means of adding the results obtained in each element.We compute where in each element we numerically propagate the initial condition at the starting point by means of the dynamical equations, in order to obtain the velocity function v 2 inside the kth element.Then a Simpson quadrature is employed to compute the integral.
To get an estimation of the error inside a given element, we consider the former v 2 velocities inside the element v k 2 and the velocities v k 1 obtained taking the derivative of the approximate solution given by the finite element method as well inside this element a linear function in our case .An estimator of the error inside the element is computed by means of where ν is the acceptability criteria, the threshold parameter control of the adaptive remeshing procedure that will be discussed and tuned in the examples of Section 4.
In order to compute a new mesh when it is not accepted, we use the Li and Bettess remeshing strategy see 13 .This strategy is based on the idea that the error distribution on an optimal mesh is uniform where ν is again the acceptability criteria, e k is the computed error on element k, M is the number of elements of the mesh, and the hat distinguishes the parameters of the new mesh to be generated.The strategy consists on finding the new length of the elements using the number of elements of the new mesh, M. According to Li and Bettess, if d denotes the dimension of the problem and m the maximum degree of the polynomials used in the interpolation for the approximate solutions inside an element k, then the number of elements that should have the new mesh is

3.6
Since we work with linear elements in dimension one, we have m 1 and d 1, and the recommended number of elements of our new mesh is

3.7
Once we have the estimation of the number of new elements, we can find the length of them by means of that, in our case, turns out to be 3.9

Simulations with Adaptive Remeshing
As it has been mentioned, in the following simulations we have located the formation in the vicinity L 2 in the Sun-Earth system.In particular we choose a halo orbit of 120000 km of z-amplitude.
We have considered two limiting cases.The first one involves no collision risk.It is known that in this case the optimal solution for each spacecraft is a bang-bang control and in this section we see that our methodology converges towards it.We note that this is the most critical case for our procedure, since the optimal maneuver is not a continuous function but it consists of two impulsive delta-v: one at departure and another one at the arrival position.The remaining nodal delta-v must be zero, and consequently this is a case where the computation of derivatives for J 1 is very ill conditioned.
The second limiting case corresponds to reconfigurations with collision risks.In this case the simple bang-bang controls would cause collision between the spacecraft, and the FEFF methodology solves the problem tending to low thrust solutions when the diameter of the mesh tends to zero.This is, the methodology can cope with both impulsive and smooth function controls selecting the optimal one for each case or part of the trajectory.
With the purpose of calibrating the acceptability parameter ν, in this section we present some simulations with reconfigurations in different situations involving and not involving collision risks.

A Bang-Bang Simulation Considering a Single Spacecraft
When the reconfiguration maneuver is not affected by collision risk for one or more spacecraft of the formation, these satellites follow independent trajectories i.e., the collision avoidance constraints in fact will not be active .So we can consider a formation just consisting of a single spacecraft to exemplify the procedure.
Let us consider in this example a shift of a single spacecraft.The reference frame for 2.2 is aligned with respect to the RTBP one but with origin on the nominal point of the base halo orbit which has been taken of 120000 km in the z-amplitude .When t 0, this point on the halo orbit corresponds to the "upper" position, this is, when it crosses the RTBP plane Y 0 with Z > 0. The initial condition for this example is taken 100 meters far from the base nominal halo orbit in the X direction, and the goal is to transfer it to a symmetrical position with respect to the halo orbit in 8 hours.This is to 100 meters in the opposite X direction performing a total shift of 200 meters during the transfer maneuver.For this particular case we know that the optimal solution is a bang-bang control with maneuvers of 0.69 cm/s at departure and arrival.
Our procedure starts with FEFF-DV2 minimizing the J 2 functional 3.1 and a trajectory with the delta-v profile of Figure 2. In fact, as we have discussed previously, since there are no collision risks, this optimal solution corresponds to the initial seed we provide for FEFF-DV2.Moreover, if we do not take into account the magnitude of the maneuvers for this particular example, the delta-v profile displayed in Figure 2 is the usual one we find in similar situations.As expected, since it does not correspond to the optimal solution of the J 1 functional 2.3 , it is not a bang-bang control.
Using this solution as the initial seed, we start the second part of the procedure corresponding to the iterative part in the schema of Figure 1.It involves FEFF-DV, the estimation of the error for the current mesh and the generation of the new one.In Figure 3 we can see the delta-v profile that we obtain after the iterations one and three; this last one is already very close to the bang-bang control the method converges after 5 iterations .
Of course in a real situation, and specially for small thrusters, the maximum value of delta-v may be constrained.In Figure 4 we show the delta-v profile for the same simulation but constraining its maximum value.Values of 0.4 cm/s and 0.3 cm/s have been chosen.We see how the methodology splits the optimal impulsive bang-bang control of about 0.69 cm/s in longer arcs at the end points of the trajectory.

Considerations and Calibration of the ν Parameter
Let us consider now the impact of the parameter ν in the performance of the procedure.We note that this parameter does not only appear in the acceptability criteria 3.4 , but it is also used to obtain the new mesh in 3.9 .
Intuitively one can expect that if we take a small value of ν, we could end up with a mesh with a big number of nodes, which turns into an optimization problem with a very large number of variables.If we take into account that to a mesh of 100 elements we have associated an optimization problem of 594 variables, we could end up with unsolvable problems in practice.On the other way round, if we use a big ν, we could end up accepting some meshes with big errors.In Table 1, we have a summary of the results obtained for different values of the parameter ν, the number of iterations needed to reach the bang-bang solution, and the number of elements after the first iteration of the methodology.We note that when ν is very small, convergence can fail.The case with ν equal to 0.0001 makes the optimal procedure awkward.When ν is 0.001, the final number of elements is greater than 1 that we know is our final target number although the maneuvers associated to the central nodes are very small.Moreover, when ν is big, there is no convergence: the final mesh contains more elements than expected because it passes the acceptability criteria before converging to the bang-bang control.For this bang-bang case, we can conclude that the best values for ν are inside the range 0.04, 0.06 .

A Simulation Using the TPF Formation
For this case we consider a configuration based on the Terrestrial Planet Finder TPF model see 1 .We assume that the satellites are initially contained in the local plane Z 0, with the interferometry baseline aligned on the X axis.The length of the baseline is 150 meters.We simulate the swap between two pairs of satellites in the baseline: each inner satellite changes its location with the outer satellite which is closest to it in position inner satellites are maneuvered to attain outer positions and vice versa as shown in Figure 5 .Again we consider 8 hours for the reconfiguration maneuver.The process of swapping positions is affected of collision risk for any radius of the sphere of influence, and simple bang-bang controls are no longer valid.We have considered a sphere of radius 10 meters, and the FEFF methodology obtains solutions of the form displayed in Figure 5.We note that the convergence of the methodology does not depend on the radius of the sphere.The reconfiguration cost increases with the radius.The final number of elements also increases with the radius.
A discussion for the parameter ν similar to the one in the previous example is also valid here: using a small ν, we can end up with a mesh with too many elements.For example, taking ν 0.0005, in the first iteration we have around 1000 elements and we do not have only the problem of having very small elements.The optimization problem has 29970 variables something that it is not desirable at all.Also, if we take a big ν, we end up with a mesh with very few elements and a big truncation error.
In Table 2 we display a summary of the results obtained for different values of the parameter ν including the number of iterations till the methodology converges Iter , and the number of elements in the first and last iterations, M 1 and M F .Due to symmetry reasons, the number of elements in each spacecraft trajectory is the same.We note that now the best values are inside the range 0.005, 0.05 and that the value ν 0.05 is appropriated for the two  cases studied.Since the shift case and the swapping case are someway the building blocks of the reconfiguration maneuver, we could suggest that ν 0.05 is a good value for the computation of reconfiguration maneuvers by means of adaptive remeshing.

Mixed Case Simulation: 3 Spacecraft
We consider in this section a case that would demand both a bang-bang and low thrust arc in the reconfiguration maneuver.The formation consists of 3 spacecraft which are in the same plane as shown in Figure 6.The reconfiguration is the swap of two spacecraft and the shift of the third one.If we perform the transfers sequentially in time, the maneuver decouples into two independent problems like the ones considered in the previous examples bang-bang plus swap .However, we are going to consider all the transfers in parallel in the same time interval, this is, with a collision risk of the three satellites in the center of the formation.
Again, we have a similar discussion about the parameter ν and the simulation results are summarized in Table 3.In this case, the values of ν that are good for our purposes are inside the range 0.002, 0.05 .

Summary of Considerations about the Value of ν
In the previous sections we have seen that a desirable value for the adaptive remeshing control parameter ν should be in the range 0.005, 0.05 .This range already gives us an idea of the value of ν that we should choose when using adaptive remeshing for the computation of reconfiguration maneuvers by means of FEFF.
We have applied the reconfiguration procedure to a test bench of 25 reconfigurations which include swaps between spacecraft located at opposite vertices of polygons 6 cases , swaps in the TPF formation 9 cases , and parallel shifts 10 cases .Different sizes and number of spacecraft, from 3 to 10, have been considered.Ten of the reconfigurations would be converging to a bang-bang solution while the other 15 would converge to low-thrust arcs when the diameter of the mesh tends to zero.We have considered our methodology taking different values of ν, and we have computed the mean of the number of iterations of the adaptive process necessary to converge.The obtained results are summarized in Table 4 and point again to the value of ν 0.05 as a convenient parameter for this kind of proximity maneuver computations.

Conclusions
This paper presents an adaptive remeshing strategy applied to a methodology to find trajectories for reconfigurations of spacecraft formations.The procedure adapts the mesh in a systematic and optimal way, and a suitable value for the parameter controlling the procedure has been found.Moreover, the methodology is robust in all the ranges of possible reconfiguration cases: from the ones that should result in a bang-bang control to the ones that should be performed with low thrust arcs.where Ω x, y, z x 2 y 2 /2 1 − μ /r 1 μ/r 2 1 − μ μ/2, μ is the mass of the small primary, and r 1 and r 2 are the distances from the spacecraft to the big and small primaries respectively.
Writing A.1 as a system of first order equations, ẋ f x , we have that x x 1 , x 2 , x 3 , x 4 , x 5 , x 6 is the state vector x, y, z, ẋ, ẏ, ż , and f is given by A.2 The reference system we consider in this paper for 2.1 is parallel to the one of the RTBP but with origin in a halo orbit; that is, it is a time-dependent translation of the synodic RTBP one given by X t x t − x h t , A.3 where x h t is the current state on the chosen halo orbit.Linearizing then ẋ f x about the halo orbit, we have Ẋ ẋh f x h Df x h X, and since x h t is a solution of ẋ f x , we obtain Ẋ Df x h X, which defines A t Df x h t in 2.1 .

, 3 . 3 where 1 e 2 2• • • e 2 M
t k and t k 1 are the ends of the kth element.From these values we have an estimation of the error of the mesh: ||e|| e 2 , and we accept the mesh when e ≤ ν u , 3.4

Figure 2 :
Figure 2: Delta-v obtained in the minimization of the J 2 functional 3.1 in a case of no collision risk.

Figure 3 :Figure 4 :
Figure 3: Delta-v profile obtained with the minimization of the J 1 functional 2.3 in a case without collision risk.a we have the result after the first iteration and b we show the result after three iterations.

Figure 5 :
Figure 5: Example of reconfiguration with collision risk: the switch of two pairs of spacecraft of the TPF formation.a shows the schema of the reconfiguration and b is a solution of the reconfiguration obtained with the FEFF methodology.

Figure 6 :
Figure 6: Example of reconfiguration with mixed bang-bang and collision risk: the swap of a pair of spacecraft and a shift.
Using this reference frame, the equations of motion of the RTBP are ẍ

Table 1 :
Number of iterations necessary to obtain the bang-bang solution of the first example depending on the parameter ν.We have indicated by "fail" the cases where the procedure does not converge and in M 1 the number of elements in the first iteration.

Table 2 :
Number of iterations and elements involved in the swapping example of TPF depending on the parameter ν.M 1 refers to the number of elements in the first iteration and M F to the final ones.

Table 3 :
Number of iterations needed to obtain the mixed solution depending on ν.

Table 4 :
Mean value of the number of iterations as a function of ν for the 25 test bench reconfigurations considered in the simulation.