Differential Drag Control Scheme for Large Constellation of Planet Satellites and On-Orbit Results
Cyrus Foster, James Mason, Vivek Vittaldev, Lawrence Leung, Vincent Beukelaers, Leon Stepan, Rob Zimmerman
IIWSCFF 17-33DIFFERENTIAL DRAG CONTROL SCHEMEFOR LARGE CONSTELLATION OF PLANET SATELLITESAND ON-ORBIT RESULTS
Cyrus Foster ∗ , James Mason † , Vivek Vittaldev ‡ , Lawrence Leung § ,Vincent Beukelaers ¶ , Leon Stepan (cid:107) and Rob Zimmerman ∗∗ A methodology is presented for the differential drag control of a large fleet ofpropulsion-less satellites deployed in the same orbit. The controller places satel-lites into a constellation with specified angular offsets and zero-relative speed.Time optimal phasing is achieved by first determining an appropriate relativeplacement, i.e. the order of the satellites. A second optimization problem is thensolved as a large coupled system to find the drag command profile required foreach satellite. The control authority is the available ratio of low-drag to high-dragballistic coefficients of the satellites when operating in their background mode.The controller is able to successfully phase constellations with up to 100 satellitesin simulations. On-orbit performance of the controller is demonstrated by phasingthe Planet Flock 2p constellation of twelve cubesats launched in June 2016 into a510 km sun-synchronous orbit.
INTRODUCTION
Due to easier access to space in recent years, there have been an increasing number of satelliteslaunched. In addition to more launches, multiple satellites are deployed into orbit in one launch.The resent launch by the Indian Space Research Organisation (ISRO), for example, carried 104satellites into orbit in February 2017. Other launches with a large number of satellites were aRussian Dnepr launch with 37 satellites, and an Orbital Antares rocket with 34 satellites onboard in2004. Such launches are optimal for large constellations of small satellites that are in the same orbitbut spread out in mean anomaly. The customized or randomized deployment velocities providedby the launcher can help in the spreading of the satellites. It is, unfortunately, not possible to firstspread the satellites and then maintain relative phasing by only using the deployment velocities. Torealize a fully spread and operationally useful constellation as shown in Figure 1, some propulsivecapability is required.A constellation launched together typically consists of small homogeneous satellites. Addingpropulsion to these satellites is an undertaking in terms of engineering, cost, and regulations. Itmight not be possible to fit the required propulsion subsystem for the desired ∆ v due to size and ∗ Orbit Mechanic, AIAA member, [email protected] † VP Missions, AIAA member, [email protected] ‡ Missions Operations Analyst, AIAA member, [email protected] § ADCS Engineer, AIAA member, [email protected] ¶ ADCS Engineer, AIAA member, [email protected] (cid:107)
Missions Operations Analyst, AIAA member, [email protected] ∗∗ Systems Engineer, AIAA member, [email protected] a r X i v : . [ phy s i c s . s p ace - ph ] J un igure 1 : A line scanner configuration of satellites uniformly spread over the orbitweight constrains. Atmospheric drag, however, can be used to provide the relative velocities re-quired to phase the constellation without additional mass if the altitude is sufficiently low. Differen-tial drag changes the drag characteristics of space objects in order to control relative accelerations.The effective surface area perpendicular to the velocity vector acts as the actuator by changing theBallistic Coefficient (BC).Differential drag was first proposed for station keeping of a pair of satellites in 1989 by usingdrag plates. The control freedom is the angle of attack of the drag plates at ◦ or ◦ and the controllaw is based on the linearized Clohessy-Wiltshire equations of motion . A similar strategy usinga drag plate and including J perturbation for satellite rendezvous is demonstrated by Bevilacqua using the Schweighart-Sedwick linearized equations of motion . Rendezvous using a combinationof differential drag for the first phase and low-thrust propulsion for the second precise phase isshown . Differential drag has also been combined with other perturbations such as Solar RadiationPressure (SRP) for formation maneuvers .Various types of controllers have been used to create trajectories such as an optimal control ap-proach , Lyapunov based controllers , and sliding mode control . The extension of the two satelliteproblem into a multiple satellite problem has been carried out previously .Differential drag has also been used for collision avoidance . For satellites using differentialdrag, uncertainties in the atmosphere models have a noticeable effect on the uncertainty propagationof the satellites . Another possible application of differential drag is to minimize the post-deployment probability of collision of nanosatellties .To date, the Aerospace Corporation has demonstrated the use of differential drag on-orbit withsmall satellites, most recently with the AeroCube-4 mission in 2013 . This three-satellite missionwas deployed into a 480 x 780 km orbit and operators were able to control along-track spacingof the satellites by selectively retracting solar panel wings. Orbcomm, with its constellation of 35satellites in a 720 km altitude circular orbit, has also used differential drag to save propellant on itsthrust-based station-keeping system .The work presented in this paper shows an end to end implementation of a differential dragcontroller that phases out and performs station keeping on a constellation of satellites. The proposedcontroller converts the phasing problem into two separate optimization problems. The phase called Reslotting is only relevant when the constellation has more than two spacecraft where the relativephasing of the satellites is determined. The atmospheric density is allowed to be time varying, as2s the case for real missions. Almost all prior work in differential drag, except for Guglielmo etal. , assumes a constant density. Intra-constellation coupling is also included when computing thecontrol. Coupling is considered by Harris and Ac¸ikmes¸e for up to five spacecraft with constantatmospheric drag. It is, however, unclear how this method scales up to 100 satellites. Li andMason treated a degenerate case of the multi-spacecraft problem limited to five satellites wherethe reference satellite performs no high drag maneuver. In addition to simulation results, on-orbitresults of the successful phasing of 12 cubesats and the ongoing phasing of 88 cubesats is shown. DIFFERENTIAL-DRAG CONTROLLER
As illustrated in Figure 2, the controller consists of three components. It first uses Cartesianposition-velocity state vectors generated by orbit determination to estimate the relative motion ofthe satellites in the along-track direction. Two optimization problems aiming to minimize the totalphasing time are then solved sequentially:1. Slot Allocator: allocating each satellite to specific slots2. Command Generator: finding appropriate high-drag windows required to reach the slots
Differential Drag Controller
Desired slot configurationStatevectors
Relative State Estimator Slot Allocator Command Generator
High-drag windows
Figure 2 : Controller produces high-drag windows from state vectors
Relative State Estimator
Since differential drag only controls the along-track dynamics of the satellite, the independent setof cartesian state vectors are reduced to relative Earth-centered angles shown in Figure 3 and theirtime derivatives.
Initial Condition
Relative angles θ are calculated between each satellite and a reference satellite.This reference can be chosen to be any of the satellites but by convention the one that starts in thelowest orbit is used. θ osc is defined as the osculating Earth-centered relative angle between each satellite i and thereference’s position vectors R , defined over [0 .. π ) in the along-track direction (Equation 1). Thehalf-plane ambiguity is solved with the reference pole vector (cid:126)N = || (cid:126)R ref × (cid:126)V ref || . It is noted thatthe reference satellite’s θ and ˙ θ are always 0 by definition. θ osc,i = angle [0 .. π ) ( (cid:126)R ref , (cid:126)R i ) (1)3 eference satellite (cid:7520) = 120° (cid:7520) = 30° (cid:7520) = -20° Figure 3 : Relative state θ reduces complexity of control problem to along-track spaceThe initial state consisting of θ and its time derivative, the angular velocity ˙ θ , are calculatedas mean elements by fitting a line through the time history of θ osc via least-squares as shown inFigure 4a. A 1-day fit period proves to be sufficient for capturing the motion of a typical LEO SSOorbit (400-600 km). To produce this time history, orbits are propagated under a 10x10 gravity andNRL-MSIS00 atmosphere force model. Control authority
To model the dynamics of satellites in high-drag, ¨ θ is also calculated, therelative angular acceleration of one satellite in high-drag mode relative to the reference satellite inlow-drag. ¨ θ is a measure of available control authority, and is calculated across the simulation timehorizon to effectively model the time-varying nature of control authority due to osculating motionof perigee and apogee and forecasted solar flux variations. This time-varying nature is importantto capture when planning months ahead, and especially so for elliptical orbits where perigee andapogee vary cyclically over month timescales due the Earth’s gravity field (Figure 4b).The mean dynamic pressure encountered by the reference satellite q ref throughout a discretizedtime period ∆ t at time k is first obtained using the atmosphere model with the latest space weatherdata. q ref ,k = 12 ρv (2)Next, the ∆ V during the discretized time period is calculated using the observed ballistic coeffi-cients BC = mC d A in low and high-drag modes, BC L and BC H respectively: ∆ V = q ref ,k (cid:16) BC L − BC H (cid:17) ∆ t (3)The mean relative velocity ˙ y between the satellites is then evaluated, using the secular along-trackcomponent of the solution to Hill’s equations for relative motion : ˙ y ≈ − y = − V = 3 q ref ,k (cid:16) BC H − BC L (cid:17) ∆ t (4)4 θ [ d e g s ] osculatingfit (a) Relative state θ and ˙ θ from osculating motion a v e r a g e a t m o s p h e r e d e n s i t y ρ [ k g / m ]
1e 13 (b) Mean atmosphere density encountered by satellite,and hence control authority for differential drag, can fluc-tuate with time
Figure 4 : Relative motion is estimated from propagated satellite state vectorsFinally, the angular acceleration ¨ θ k is evaluated during this discretized period, converting to an-gular coordinates using the reference satellite’s mean semi-major axis at epoch a ref : ¨ θ k = ˙ y ∆ ta ref = 3 q ref ,k a ref (cid:16) BC H − BC L (cid:17) (5) Dynamics
The along-track dynamics of a constellation of n satellites subject to time-discretizedbinary control inputs (low-drag or high-drag) can now be summarized by n − systems of Equation6 for each satellite-reference pair. Note these systems are coupled since they all share the reference’scommands u ref . (cid:20) θ ˙ θ (cid:21) k +1 = (cid:20) t (cid:21) (cid:20) θ ˙ θ (cid:21) k + (cid:20) B k ( u k ) (cid:21) (6a)where the control inputs B k are evaluated over the discretization period ∆ t with the combinedangular acceleration from the two satellites: B k ( u k ) = ( − ¨ θ ref,k u ref,k + ¨ θ sat,k u sat,k )∆ t (6b)and high-drag commands u i,k are defined by: u i,k = (cid:40) if i is in low-drag at t k if i is in high-drag at t k (6c)5 Θ ( d e g s ) phasing0 10 20 30 40 50 60 70 80 90Days (a) random slot order: ∆ t phase = 90 days Θ ( d e g s ) phasing0 10 20 30 40 50 60 70 80 90Days (b) optimal slot order: ∆ t phase = 70 days Figure 5 : Allocating slots judiciously significantly lowers the time required to phase a constellation.
Slot Allocator
With the control problem reduced to optimization in relative-motion θ space, satellite orderingrelative to one another is first determined to minimize predicted phasing time. The Dove satellitesare interchangeable, so optimization is targeted to seek the lowest phasing time that achieves thedesired slot vector irrespective of satellite order. As shown in Figure 5, the judicious assignmentof satellites to slots is important for minimizing the time to phase of the constellation. Figure 5ashows a phasing attempt when target slots are assigned randomly, while Figure 5b starts with thesame initial condition but slots have been allocated optimally to minimize time to phase. Slots
A vector of slots is defined, where a slot is a desired θ f for the final state of slot with ˙ θ f = 0 . Various examples of slot configurations are shown in Figure 6. Equally-spaced satellitesmake sense for minimizing imaging swath overlap and ground station conflicts, while the fixed-spaced configuration is ideal for realizing a line scanner with adjacent swaths. Fixed-spaced e.g. 25° spacing
Custom e.g. [0°, 15°, 180°, 195°][0°, 25°, 50°, 75°] [0°, 15°, 180°, 195°]
Resultant slot vector: Equally-spaced [0°, 90°, 180°, 270°]
Figure 6 : Slot vector describes desired relative placement of satellites6 lip-flop solution
Since the objective is to minimize phasing time, there is a need to estimate theamount of time required to achieve a certain slot configuration (cid:2) θ ˙ θ (cid:3) Tf given the initial condition (cid:2) θ ˙ θ (cid:3) T of each non-reference satellite. A method is introduced for solving the phasing time of atwo-satellite constellation open-loop; to be used subsequently in the slot allocation optimizer.The time-optimal solution for phasing two-satellites always has two phases (Figure 7): one satel-lite high-drags for a duration ∆ t A while the other low-drags (flip), then the roles reverse for anotherduration ∆ t B (flop). The solution to this control problem (duration of each phase, which satellitestarts first) can be found analytically if assuming time-invariant ¨ θ (see Appendix B). The generalproblem can also be solved numerically via a line search as outlined in Algorithm 1, where therelative motion is propagated with the discretized equations of motion from Equation 6 to accountfor time-varying ¨ θ . Θ ( d e g s ) Referencehigh-drags Satellitehigh-drags low-draghigh-drag
Figure 7 : Flip-flop solution provides the time-optimal schedule for phasing two satellites
Algorithm 1
Flip-flop solution procedure S OLVE TWO - SATELLITE FLIP - FLOP PROBLEM try both order combinations for first sat to high-drag (one has a bracketed root in line search) while line searching to find ∆ t A such that θ f = θ f,desired do prop via Eq 6 first sat in high-drag for ∆ t A then second in high-drag until ˙ θ f = ˙ θ f,desired end while output ∆ t phase = ∆ t A + ∆ t B , and order in which satellites high-drag. end procedure Optimizer
The optimal slot configuration is solved for by selecting the order which minimizesthe largest ∆ t phase of each reference-satellite pair. Algorithm 2 describes an implementation usingsimulated annealing wherein random perturbations are made to the slot configuration and kept ifit satisfies an acceptance probability described in equation 7. It is found that values of k max = 10 and t = 100 produce good results for typical differential drag scenarios.7 Θ ( d e g s ) low-draghigh-drag (a) without commands Θ ( d e g s ) phasing station-keeping low-draghigh-drag (b) with commands Figure 8 : Time-discretized high-drag commands are assigned to achieve desired slots
Algorithm 2
Slot Allocator procedure O BTAIN SLOT CONFIGURATION FOR MINIMUM TIME PHASING for k through k max do create slot new from slot by randomly swapping two satellites compute ∆ t to phase for all reference-sat pairs with Algorithm 1 if P (∆ t max,old , ∆ t max,new , k, k max ) ≥ random (0 , then slot ← slot new end if end for end procedure P ( cost , cost new , k, k max ) = (cid:40) if cost new < cost exp ( cost − cost new t ) if cost new ≥ cost with t = t (1 − kk max ) (7)If the condition ∆ t max,old = ∆ t max,new is reached, because the reference-satellite pair with themaximum phasing time is the same in the nominal and perturbed slot configurations, one thenbacks up to comparing the second greatest ∆ t phase (or next greatest if still equal).When computing ∆ t to phase to a particular slot, one also solves for which π offset to use bytrying several neighboring π multiples and selecting the one with lowest phasing time. Command Generator
Now that a final desired state θ f is determined for each satellite, a set of high-drag commands isneeded to guide the satellites to their desired slots in minimum time (Figure 8).High-drag commands are discretized for simulation with Equation 6 across a time horizon thatsufficiently captures the required settling time as estimated by the flip-flop solutions. The initial8uess for the command matrix u comes from the superposition of commands from flip-flop solutionsto each reference-satellite sub-problem, or alternatively from the commands of the last time thecontroller was run. u is then perturbed to converge satellites to their assigned slots via simulatedannealing per the procedure described in Algorithm 3. Algorithm 3
Command generator procedure O BTAIN HIGH - DRAG COMMANDS TO TARGET SLOTS IN MINIMUM TIME for k through k max do create u new with randomly flipped command u sat = i,t = k compute finals states θ i,f under u new with eq. 6 compute cost new = (cid:80) i ( θ i,f − θ i,desired ) if P ( cost old , cost new , k, k max ) ≥ random (0 , then u ← u new end if end for end procedure During the station-keeping phase, the default high-drag mode is replaced with one that performshigh-drag only a fraction of the time (e.g. 25%) to minimize the lifetime lost to drag since fullcontrol authority is no longer required. This also has the desireable side effect of shortening thelimit cycle (station-keeping phasing error) when the time discretization period is kept constant.
RESULTS
Results from the simulation and on-orbit implementation of the controller presented in Sec-tion are presented here. All simulations and results use Planet’s Dove satellites. The variousattitude modes of Planet’s Dove satellites exhibit a large difference in cross-sectional area thanks totheir deployed solar panels, as shown in Figure 9.Although there is a 10:1 ratio between solar-panel and telescope facing cross-sectional areas, onlya 3:1 ratio in orbit-derived ballistic coefficients (BC) is observed in reality. This reduction in controlauthority is due to several factors: • Satellite duty-cycle: a satellite only executes the desired low-drag or high-drag mode when itis neither imaging nor downlinking. In imaging mode, the satellite presents a cross-sectionalarea between low and high-drag areas, while during downlinks the satellite tracks the groundstation with its antenna boresight (co-aligned with telescope). • Attitude pointing accuracy: attitude errors effectively reduce the ratio of high to low-dragcross-section areas. For example when not imaging or downlinking, the satellites are notusing the star camera for attitude determination, resulting in deviations from modeled area. • Skin friction effects: even in a perfect low-drag attitude, the large solar panel faces parallel tothe velocity vector still contribute to drag via skin friction effects , resulting in the low-dragmode having a higher drag coefficient c d . 9 ow-drag modecamera facing200 cm high-drag modesolar panel facing1950 cm imaging modeside-panel facing370 cm (a) Orthographic projections with cross-sectional areas. (b) High-drag and low-drag attitudes Figure 9 : Attitude modes of Dove satellite enable large drag area ratios
Performance of Controller
The main contributions of the proposed controller is the slot allocator and solving for the trajec-tories of all satellites as one highly coupled system. The coupling occurs due to the various flip-flopsolutions possible requiring conflicting high-drag windows for the reference satellite.Although the current controller reliably produces a robust solution, there are a few possible forimprovements: a one-step optimizer and second-order optimizations. A source of future work is tosolve the two optimization steps in one, by both assigning slots and generating commands in a singleoptimization loop. This would account for inter satellite-pair coupling effects, but it is currentlycomputationally prohibitive to nest the command generator in the slot allocator. Finally, it wouldbe of interest to implement second-order objective functions such as minimizing clustering duringearly stages of phasing, and minimizing station-keeping control effort given a specified control box.
In-Space Results
For satellite operations, all estimation and control for differential drag is performed on the ground(Figure 10). GPS data is regularly downloaded during X-band passes and ground-based orbit deter-mination maintains state vector ephemerides for each satellite. The differential drag controller thenuses these state vectors and the desired slot configuration to produce a set of high-drag commandsthat are uploaded to the satellite. For Planet’s operational orbits in the neighborhood of 500 km alti-tude, it is found that a controller update and time discretization of 1-day are sufficient for achievingdesired performance.The differential drag controller described in this paper was applied to Flock 2p, a set of 12 satel-lites launched into a circular 510 km sun-synchronous orbit in June 2016. The satellites were de-ployed with a non-uniform total spread of 0.5 m/s along-track ∆ V (deployer ejection speed is 1 m/s,along-track spread results from upper stage attitude schedule), and the controller allocated high-dragcommands that initially increased this spread before converging to the desired slot configuration.The orbit-derived relative motion and commanded high-drag windows are shown in Figure 11. The10 rbitDetermination Differential Drag Controller State vectorsDesired slot configurationGPS data High-drag windowsSpace segmentGround segment
Satellites
Figure 10 : Orbit determination and differential drag controller runs on the ground - - - - - - - - - - - - - - - - - - - - - - R e l a t i v e a n g l e ( d e g s ) low-draghigh-drag Figure 11 : On-orbit results from Flock 2psame information is presented in a polar plot in Figure 14 of Appendix A. This relative motion plotcan be recreated using JSpOC TLEs (CATIDs 41606, 41608-41618) or publicly-available Planet-produced ephemerides ∗ . The plan for Flock 2p was to achieve an equally-spaced constellation tominimize swath overlap and antenna conflicts, except for a pair of satellites that would maintain ∗ http://ephemerides.planet-labs.com/ ug 2016 Sep 2016 Oct 2016 Nov 2016 Dec 2016020406080100 B a lli s t i c c o e ff i c i e n t ( k g / m ) actual low-dragexpected low-dragactual high-dragexpected high-drag (a) Ballistic Coefficients (BC) Aug 2016 Sep 2016 Oct 2016 Nov 2016 Dec 20160.000.010.020.030.040.050.060.07 C o n t r o l a u t h o r i t y ( d e g r ee s / d a y ) actualexpected (b) Control authority Figure 12 : Expected vs. actual atmosphere-dependend performance . ◦ spacing to demonstrate the line scanner imaging strategy.The primary lesson learned from automated on-orbit control of Flock 2p is to account for thesignificant unmodeled variations in effective BC. These fluctuations have been treated extensivelyby others , and its specifics as encountered by Planet’s differential drag phasing are discussedbelow. Figure 12a shows the expected versus actual BCs; the expected values were pre-launchestimates based on an older Flock launched into a 600 km orbit in 2014. Figure 12b shows theresulting available control authority with pre-launch BCs using the solar flux prediction data at thetime, and the actual control authority available from post-processed orbit data. One can see that theactual control authority available was significantly less than that predicted, by up to 50% 6 monthsafter launch. This underestimation of control authority is responsible for the overshoot observed inFigure 11; and caused some churn with the slot allocator as it attempted to continuously target slotsusing its optimistic estimate of control authority.Even though the satellites’ attitude control systems were consistent at maintaining low and high-drag modes, the resulting BCs exhibit large periodic fluctuations. The same fluctuating behaviorhas also been observed across year-long time-scales on satellites that are tumbling at rates muchshorter than orbit fit-spans (therefore effectively presenting a constant cross-sectional area). Onewould expect the BCs in various consistent modes to be time-invariant when using a-posteriorisolar flux data and industry standard atmosphere models like MSIS00 or Jacchia-Bowman 2008 ,so the fluctuating behavior of orbit-derived BCs is attributed to unmodeled atmospheric densityvariations due to solar phenomena (the period of the fluctuations also matches the Sun’s synodicperiod of 26 days). Atmosphere models have evolved in complexity, but they still fundamentallyrely on a handful or less of channels of solar data. Since the Sun is the greatest contributor tothese atmospheric density fluctuations, there must be potential for improved atmosphere modelsthat make use of more data from the Sun, perhaps 10’s or 100’s channels and at higher frequenciesfrom various measurement platforms that were not historically available.The mitigation was to continuously update the controller’s estimate of low and high-drag BCsbased on observed control authority from satellites in those modes.12 uture Contellations Planet recently launched Flock 3p on February 14th 2017 into a 505 km altitude orbit. It consistsof 88 Doves that were individually deployed and targeted to achieve a roughly uniform distributionof along-track ∆ V ’s. The total initial spread is 2 m/s, spanning two groups of satellites, with agap in the middle consisting of 13 non-Planet satellites. Figure 13 shows the simulated behaviorof Flock 3p from deployment through station-keeping. The same information is presented in apolar plot in Figure 15 of Appendix A. The simulation factors in the fact that not all satellites areeligible to perform differential drag maneuvers right-away, as commissioning of specific satellitesare staggered over the first 80 days. - - - - - - - - - - - - - - - - - - - - - - R e l a t i v e a n g l e ( d e g s ) low-draghigh-drag Figure 13 : Simulation of Flock 3p phasing and station-keeping
CONCLUSION
This paper details a controller design for phasing and station-keeping large fleets of satellites withonly differential drag. The controller works in three steps: relative motion estimation, slot allocatingand high-drag command generation for the entire couple system.Relative motion is first estimated in the along-track direction to obtain the mean angular separa-tion angle and speed between each satellite and a designated reference. The satellites are assumedto be commandable into discrete low and high-drag modes, and angular acceleration for each modeis evaluated across a time horizon of interest to capture variations in mean atmospheric density.Two optimization problems are then solved sequentially: the first to assign each satellite to a spe-cific constellation slot and the second to generate time-discretized commands to guide each satellite13o its slot using a coupled system. Both optimization problems seek to minimize the required phas-ing time given the available control authority, and are solved with simulated annealing, a methodthat lends well to the discretized nature of the problem. They also make use of the solution to thetwo-satellite sub-problem to either estimate the required phasing time or provide an initial guess tothe general multi-satellite problem. Although both components are required for an optimal solution,the coupling of the system has a larger influence on the optimality and feasibility of the solutions.The closed-loop performance of the controller is demonstrated with an operational fleet of twelvesatellites, and it is also being applied to another fleet of 88 satellites that was recently launched inFebruary 2017.
ACKNOWLEDGMENTS
The authors would like to acknowledge Planet’s Missions and Spacecraft Design teams for puttingtogether all the pieces needed to execute the automated differential drag operations described in thepaper.
APPENDIXPolar Plots
Figures 14 and 15 are polar projections of Figures 11 and 13 respectively. Radial axis shows time(in days), while relative angle is on the θ axis. These polar projections better illustrate the relativeangles that wrap around, at the cost of a time axis that is harder to read. low-draghigh-drag Figure 14 : On-orbit results from Flock 2p14 ow-draghigh-drag
Figure 15 : Simulation of Flock 3p phasing and station-keeping
Analytic solution for two satellite case
This section presents the analytical solution to the time-optimal open-loop differential drag con-trol of two satellites under time-invariant control authority ¨ θ . The initial relative motion is given by (cid:2) θ ˙ θ (cid:3) T and the final condition by (cid:2) θ ˙ θ (cid:3) Tf .The solution consists of two phases (Figure 7): one satellite high-drags for a duration ∆ t A whilethe other low-drags, then the roles reverse for another duration ∆ t B . If the reference high-dragsfirst, the effective control authority during phases A and B are ¨ θ A = − ¨ θ and ¨ θ B = ¨ θ respectively.Similarly if the satellite high-drags first, then one has ¨ θ A = ¨ θ and ¨ θ B = − ¨ θ respectively.The relative motion at the mid-point between phases A and B is expressed as (cid:2) θ ˙ θ (cid:3) Tm : θ m = θ + ˙ θ ∆ t A + 12 ¨ θ A ∆ t A ˙ θ m = ˙ θ + ¨ θ A ∆ t A (8)The final condition after phase B is then expressed as: θ f = θ m + ˙ θ m ∆ t B + 12 ¨ θ B ∆ t B ˙ θ f = ˙ θ m + ¨ θ B ∆ t B (9) ∆ t A and ∆ t B are solved for from equations 8 and 9 by first eliminating θ m and ˙ θ m , and intro-ducing the short hands ∆ θ = θ f − θ and ∆ ˙ θ = ˙ θ f − ˙ θ , to reveal ∆ t B as a function of ∆ t A :15 t B = ∆ ˙ θ − ¨ θ A ∆ t A ¨ θ B (10)and ∆ t A from the quadratic: (cid:34)
12 ¨ θ A (cid:18) ¨ θ A ¨ θ B − (cid:19)(cid:35) ∆ t A + (cid:34) ˙ θ (cid:18) ¨ θ A ¨ θ B − (cid:19)(cid:35) ∆ t A + (cid:34) ∆ θ − ∆ ˙ θ ¨ θ B (cid:18) ˙ θ + 12 ∆ ˙ θ (cid:19)(cid:35) = 0 (11)There are four possible solutions to the ( ∆ t A , ∆ t B ) couple given the ambiguity of which satelliteshould high-drag first, and the two roots to Equation 11, but only one couple has positive real values,the physical solution to the problem. REFERENCES [1] C. Leonard, W. Hollister, and E. Bergmann, “Orbital Formationkeeping with DifferentialDrag,”
AIAA Journal of Guidance, Control, and Dynamics , Vol. 12, No. 1, 1989, pp. 108–113, 10.2514/3.20374.[2] W. H. Clohessy and R. S. Wiltshire, “Terminal Guidance System for Satellite Rendezvous,”
Journal of the Aerospace Sciences , Vol. 27, No. 9, 1960, pp. 653 – 658, 10.2514/8.8704.[3] R. Bevilacqua and M. Romano, “Rendezvous Maneuvers of Multiple Spacecraft Using Dif-ferential Drag Under J2 Perturbation,”
Journal of Guidance, Control, and Dynamics , Vol. 31,No. 6, 2008, pp. 1595–1607, 10.2514/1.36362.[4] S. A. Schweighart and R. J. Sedwick, “High-Fidelity Linearized J Model for Satellite Forma-tion Flight,”
Journal of Guidance, Control, and Dynamics , Vol. 25, No. 6, 2002, pp. 1073–1080, 10.2514/2.4986.[5] R. Bevilacqua, J. S. Hall, and M. Romano, “Multiple spacecraft rendezvous maneuvers bydifferential drag and low thrust engines,”
Celestial Mechanics and Dynamical Astronomy ,Vol. 106, No. 1, 2009, p. 69, 10.1007/s10569-009-9240-3.[6] D. Spiller, F. Curti, and C. Circi, “Minimum-Time Reconfiguration Maneuvers of SatelliteFormations Using Perturbation Forces,”
Journal of Guidance, Control, and Dynamics , Aheadof print, 2017, 10.2514/1.G002382.[7] L. Dell’Elce and G. Kerschen, “Optimal propellantless rendez-vous using differential drag,”
Acta Astronautica , Vol. 109, 2015, pp. 112 – 123, 10.1016/j.actaastro.2015.01.011.[8] D. P´erez and R. Bevilacqua, “Lyapunov-Based Adaptive Feedback for Spacecraft PlanarRelative Maneuvering via Differential Drag,”
Journal of Guidance, Control, and Dynamics ,Vol. 37, No. 5, 2014, pp. 1678–1684, 10.2514/1.G000191.[9] S. Varma and K. D. Kumar, “Multiple Satellite Formation Flying Using Differential Aero-dynamic Drag,”
Journal of Spacecraft and Rockets , Vol. 49, No. 2, 2012, pp. 325–336,10.2514/1.52395. 1610] M. W. Harris and B. Ac¸ikmes¸e, “Minimum time rendezvous of multiple spacecraft using dif-ferential drag,”
Journal of Guidance, Control, and Dynamics , Vol. 37, 3 2014, pp. 365–373,10.2514/1.61505.[11] B. S. Kumar, A. Ng, K. Yoshihara, and A. D. Ruiter, “Differential Drag as a Means ofSpacecraft Formation Control,” , March 2007, pp. 1–9,10.1109/AERO.2007.352790.[12] M. Horsley, “An investigation into using differential drag for controlling a formation of Cube-Sats,”
Advanced Maui Optical and Space Surveillance Technologies Conference , Sept. 2011,p. E28.[13] D. Mishne and E. Edlerman, “Collision-Avoidance Maneuver of Satellites Using Drag andSolar Radiation Pressure,”
Journal of Guidance, Control, and Dynamics , Ahead of print, 2017,10.2514/1.G002376.[14] X. Huang, Y. Yan, and Y. Zhou, “Underactuated spacecraft formation reconfigura-tion with collision avoidance,”
Acta Astronautica , Vol. 131, 2017, pp. 166 – 181,http://doi.org/10.1016/j.actaastro.2016.11.037.[15] L. Mazal, D. P´erez, R. Bevilacqua, and C. F., “Spacecraft Rendezvous by Differential DragUnder Uncertainties,”
Journal of Guidance, Control, and Dynamics , Vol. 39, August 2016,pp. 1721–1733, 10.2514/1.G001785.[16] O. Ben-Yaacov, A. Ivantsov, and P. Gurfil, “Covariance analysis of differential drag-based satellite cluster flight,”
Acta Astronautica , Vol. 123, 2016, pp. 387 – 396,10.1016/j.actaastro.2015.12.035.[17] J. A. Atchison and A. Q. Rogers, “Operational Methodology for Large-Scale Deployment ofNanosatellites into Low Earth Orbit,”
Journal of Spacecraft and Rockets , Vol. 53, No. 5, 2016,pp. 799–810, 10.2514/1.A33361.[18] J. Gangestad, B. Hardy, and D. Hinkley, “Operations, Orbit Determination, and FormationControl of the AeroCube-4 CubeSats,” ,No. SSC13-X-4, Logan, Utah, 2013.[19] T. D. Maclay and C. Tuttle, “Satellite stationkeeping of the ORBCOMM constellation via ac-tive control of atmospheric drag: operations, constraints, and performance,”
AAS/AIAA SpaceFlight Mechaics Meeting , No. AAS 05-152, Copper Mountain, CO, January 2005.[20] D. Guglielmo, D. P´erez, R. Bevilacqua, and L. Mazal, “Spacecraft relative guidance via spatio-temporal resolution in atmospheric density forecasting,”
Acta Astronautica , Vol. 129, 2016,pp. 32 – 43, 10.1016/j.actaastro.2016.08.016.[21] A. S. Li and J. Mason, “Optimal Utility of Satellite Constelllation Separation with Differen-tial Drag,”
AIAA/AAS Astrodynamics Specialist Conference, AIAA SPACE Forum , No. AIAA2014-4112, San Diego, CA, August 2014, 10.2514/6.2014-4112.[22] D. Vallado,
Fundamentals of Astrodynamics and Applications . Microcosm Press and Springer,3 ed., 2007. page 387. 1723] A. Khachaturyan, S. Semenovskaya, and B. Vainshtein, “Statistical-Thermodynamic Ap-proach to Determination of Structure Amplitude Phases,”
Soviet Physics, Crystallography ,Vol. 24, No. 5, 1979, pp. 519–524.[24] G. Koppenwallner, “Satellite Aerodynamics and Determination of Thermospheric Density andWind,”
American Institute of Physics Conference Series , Vol. 1333 of
American Institute ofPhysics Conference Series , May 2011, pp. 1307–1312, 10.1063/1.3562824.[25] D. A. Vallado and D. Finkleman, “A critical assessment of satellite drag and at-mospheric density modeling,”
Acta Astronautica , Vol. 95, 2014, pp. 141 – 165,http://doi.org/10.1016/j.actaastro.2013.10.005.[26] B. R. Bowman, W. K. Tobiska, F. Marcos, and C. Huang, “The Thermospheric Density ModelJB2008 using New EUV Solar and Geomagnetic Indices,” ,Vol. 37 of