Abstract
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.
1. 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 computation of proximity maneuvers by means of artificial potential functions in the works of Badawy and McInnes [6], rotation techniques introduced by Wang et al. [7, 8], or the FEFF methodology based on a finite element implementation introduced by Garcia-Taberner and Masdemont in [9].
In this paper we consider the FEFF methodology (Finite Elements for Formation Flight) for the reconfiguration and proximity maneuvering of a formation about a libration point. With this methodology, individual trajectories of the spacecraft inside the formation are computed by means of solving a direct optimal control problem which is formulated in terms of the discretization given by a finite element procedure. We briefly summarize this procedure in Section 2. Then in Section 3 we focus on an adaptive remeshing strategy that produces optimal discretizations for the trajectories, in the sense that the error in the trajectories obtained is kept below a given threshold with the coarsest possible mesh. Finally, the paper ends with some numerical examples and conclusions.
2. 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 and is a very convenient place for space observatories ( for the Sun and for deep space). In this paper we consider a formation of spacecraft located in a halo orbit about . We assume that the formation is made of spacecraft flying with a particular pattern and our objective is to perform a reconfiguration in a fixed time . 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 is a matrix and refers to the state of the spacecraft (see The appendix). Usually the origin of the reference frame for the coordinates is the nominal point on the base halo orbit at time 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: Here and stand for the initial and final states of the th spacecraft inside the formation, and , 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 which we consider for the reconfiguration is split in 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 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, is the maneuver applied at the th node of trajectory, and and 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 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.
3. 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 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 independent problems, where we compute the trajectory which minimizes 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 ).
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 . Let us call FEFF-DV1 the formulation that uses the FEFF methodology to optimize the functional 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 , 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 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, , 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 inside the th 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 velocities inside the element () and the velocities 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 and are the ends of the th element. From these values we have an estimation of the error of the mesh: , and we accept the mesh when 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, is the computed error on element , 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, . According to Li and Bettess, if denotes the dimension of the problem and the maximum degree of the polynomials used in the interpolation for the approximate solutions inside an element , then the number of elements that should have the new mesh is Since we work with linear elements in dimension one, we have and , and the recommended number of elements of our new mesh is 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
4. Simulations with Adaptive Remeshing
As it has been mentioned, in the following simulations we have located the formation in the vicinity in the Sun-Earth system. In particular we choose a halo orbit of 120000 km of -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 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.
4.1. 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 -amplitude). When , this point on the halo orbit corresponds to the “upper” position, this is, when it crosses the RTBP plane with . The initial condition for this example is taken 100 meters far from the base nominal halo orbit in the 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 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 functional (3.1) and obtaining 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 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).
(a)
(b)
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.
(a)
(b)
4.1.1. 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 .
4.2. 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 , with the interferometry baseline aligned on the 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)
(b)
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 , 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, and . 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 and that the value 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 is a good value for the computation of reconfiguration maneuvers by means of adaptive remeshing.
4.3. 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 .
4.4. 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 . 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 as a convenient parameter for this kind of proximity maneuver computations.
5. 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.
Appendix
Let us consider the usual restricted three-body problem (RTBP) in synodic coordinates, where the unit of mass and length is taken such as the sum of the masses of the primaries and the distance between the primaries is equal to 1, and the unit of time is taken such as the period of the primaries is equal to . In this synodic coordinate system, the origin is located on the center of mass and the axis is defined by the line of the two primaries, from the smallest primary to the larger one. The axis is normal to the rotation plane, in the direction of the angular momentum of the primaries, and the axis is chosen orthogonal to the previous ones in order to have a positively orientated coordinate system.
Using this reference frame, the equations of motion of the RTBP are where , is the mass of the small primary, and and are the distances from the spacecraft to the big and small primaries respectively.
Writing (A.1) as a system of first order equations, , we have that is the state vector , and is given by
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 where is the current state on the chosen halo orbit.
Linearizing then about the halo orbit, we have , and since is a solution of , we obtain , which defines in (2.1).
Acknowledgments
This research has been supported by the MICINN-FEDER Grant, MTM2009-06973, and the CUR-DIUE Grant 2009SGR859.