Small Satellite Constellation Separation using Linear Programming based Differential Drag Commands

Small Satellite Constellation Separation using Linear Pr ogramming based Differ ential Drag Commands Emmanuel Sin 1 , Murat Arcak 2 and Andre w Packard 3 Abstract — W e study the optimal control of an arbitrarily large constellation of small satellites operating in low Earth orbit. Simulating the lack of on-board propulsion, we limit our actuation to the use of differential drag maneuvers to make in-plane changes to the satellite orbits. W e propose an efficient method to separate a cluster of satellites into a desired constellation shape while respecting actuation constraints and maximizing the operational lifetime of the constellation. By posing the problem as a linear pr ogram, we solve for the optimal drag commands for each of the satellites on a daily basis with a shrinking-horizon model predicti ve control approach. W e then apply this control strategy in a nonlinear orbital dynamics simulation with a simple, varying atmospheric density model. W e demonstrate the ability to contr ol a cluster of 100+ satellites starting at the same initial conditions in a circular low Earth orbit to form an equally spaced constellation (with a relative angular separation error tolerance of one-tenth a degree). The constellation separation task can be executed in 71 days, a time frame that is competitiv e for the state-of-the-practice. This method allows us to trade the time r equired to con verge to the desired constellation with a sacrifice in the overall constellation lifetime, measured as the maximum altitude loss experienced by one of the satellites in the gr oup after the separation maneuvers. I . I N T RO D U C T I O N Cubesats are miniature satellites comprised of one or more 10 × 10 × 10 cm cubic units. Due to the use of commercial off-the-shelf components and a widely accepted reference de- sign, cubesats ha ve become a standard platform for research [1]. As the cost of both satellite manufacturing and launch services decrease, new commercial and scientific applications of cubesats will emerge [2]. In particular , ne w opportunities may arise from the use of cubesats in lar ge-scale, coordinated constellations [3]. Coordinated groups of small satellites can enable mission objectiv es that may be difficult or impos- sible with single, monolithic satellites (e.g., high-cadence or multipoint measurements, communication relays). In fact, small satellite constellations have already been launched for commercial purposes, such as Earth imaging [4], and plans hav e been proposed for their use in providing space-based Internet [5]-[6]. Howe ver , as the number of satellites increases, it becomes more dif ficult to operate the constellation. Satellite constel- 1 Emmanuel Sin is a Graduate Student of the Department of Mechanical Engineering at the Univ ersity of California, Berkeley emansin@berkeley.edu 2 Murat Arcak is a Professor of the Department of Electrical Engi- neering and Computer Sciences at the Univ ersity of California, Berkeley arcak@berkeley.edu 3 Andrew Packard is a Professor of the Department of Mechanical Engineering at the Univ ersity of California, Berkeley apackard@berkeley.edu lation maneuvers can be divided into three broad categories: 1) initial acquisition of desired formation, 2) station-keeping of desired orbital positions and motion in the presence of disturbances, and 3) reconfiguration to a different desired formation. The current state-of-practice requires the control of each individual spacecraft from a command center on the ground that monitors the motion of each satellite. Thus, there exists the need for a centralized and optimal method of controlling satellite constellations from the ground. Cubesats typically employ a limited actuator suite com- pared to their larger cousins and may lack propulsi ve thrusters or other actuators, resulting in limited control authority . In such cases, a cubesat may ha ve to rely on “passiv e” means to make orbit maneuvers. One well-studied method called dif ferential drag [7]-[8] employs gyroscopic actuators, such as reaction wheels, to not only change a satellite’ s orientation but also its orbital in-plane motion, such as its altitude and angular speed. By changing the satellite’ s cross-sectional area exposed to the incident air molecules in the atmosphere of low Earth orbit, a varying drag force can be applied in the opposite direction of the satellite’ s orbital velocity . In what is commonly referred to as the drag paradox [9]-[10], this drag force makes the satellite fall in altitude but also increase in angular speed. By causing a differential in the angular speeds of satellites that are in the same orbit, their relativ e positions can be changed. Ho wever , differential drag maneuvers that increase a satellite’ s velocity also result in a proportional loss in its altitude, effecti vely reducing its time in orbit. Hence, these differential drag inputs must be applied sparingly so as to not unnecessarily reduce the lifetime of the satellites. Planet is a San Francisco based Earth imaging company that has successfully used dif ferential drag to form a “line- scanner” constellation using 3U (30 × 10 × 10) cubesats. The constellation allo ws Planet to produce a complete image of the Earth’ s surface e very day . Planet describes the use of a bang-bang control approach that commands each satellite to enter either a high-drag or a low-drag windo w at a certain time and for a specific duration [11]-[12]. This control strategy is considered to be operationally simpler compared to other approaches where the optimal drag area commands may take on any v alue within the continuous range between low-drag and high-drag. Ho wever , as the performance of commercially-av ailable ADCS (attitude determination and control system) sensors and actuators for cubesats increase, it may be of value to in vestigate dif ferent control strategies. The subsequent sections are organized as follo ws: Section II be gins by introducing the dynamical models used in this paper . In Section III we formulate an optimization problem that we then transform into a linear program to solve for the optimal inputs. Section IV shows our results of applying the linear program on the nonlinear orbital dynamics. W e end with Section V where we state concluding remarks. I I . DY NA M I C A L M O D E L S In this section, we first describe the nonlinear “truth” model that we use for simulating satellite orbital dynamics. W e then introduce the approximate discrete-time model to be used in our linear program. A. Orbital Dynamical Model T o determine the motion of a satellite orbiting the Earth, we start with the two-body problem where we assume that the barycenter of the system is co-located with the center of a spherically , symmetric Earth (i.e., the mass of the satellite is ne gligible). The satellite’ s motion can be described by the following second-order ordinary differential equation [13]: ¨ ~ r = − µ E | ~ r | 3 ~ r + ~ a pert urb (1) where in the first term, ~ r is the position vector pointing from the center of the Earth to the satellite and µ E is the gravitational parameter of the Earth (i.e., gravitational constant, G , multiplied by the mass of the Earth). While the first term represents the gravitational force ex erted on the satellite by the Earth, the second term, ~ a pert urb , represents the specific forces due to perturbations. Examples of such perturbations include the gravitational effects caused by the oblateness of the Earth, gravity from third bodies (e.g., Moon, Sun, other planets), solar radiation pressure and atmospheric drag. B. Atmospheric Drag P erturbation Model W e now describe the atmospheric drag force model from which we deri ve our control authority . As a satellite mov es along its orbit, it experiences an atmospheric drag force that acts against its velocity relative to the atmosphere. The mass- specific acceleration due to this atmospheric drag force is giv en by the equation [14]: ~ a at mdr a g = − 1 2 C D A m ρ | ~ v rel | · ~ v rel (2) where: C D : satellite drag coefficient A : surface area exposed to incident stream m : satellite mass ρ : atmospheric density at the satellite position ~ v rel : velocity of satellite relative to the atmosphere For the simulations in this paper , we use the same constant values used by Li and Mason [11] for the satellite’ s drag coefficient, mass, maximum and minimum drag surface areas (i.e., we use ¯ C D , ¯ m , and A min ≤ A ≤ A max ). For atmospheric density , we use a simplified v ersion of the Harris-Priester model, as described by Montenbruck [14], that does not consider diurnal effects. The relativ e velocity of the satellite with respect to the atmosphere is approximated based on the assumption that the atmosphere rotates with the same velocity as that of the Earth’ s rotation [14]: ~ v rel = ~ v sat − ~ ω E × ~ r sat (3) where ~ v sat is the satellite velocity vector , ~ r sat is the satellite position vector , and ~ w E is the Earth’ s angular velocity about its axis. C. Simulation Model In this subsection we dev elop a simple orbital dynamical model for simulation purposes. Since two-body motion is planar in an Earth-centered inertial frame, we begin by using polar coordinates to represent the satellite orbital kinematics in the plane. ~ r = r e r (4a) ˙ ~ r = ˙ re r + r ˙ θ e θ (4b) ¨ ~ r =  ¨ r − r ˙ θ 2  e r +  2 ˙ r ˙ θ + r ¨ θ  e θ (4c) W e denote the magnitude of the radial position with r and the angular position with θ . W e use e r , e θ , and e N as the unit vectors in the radial, tangential, and normal directions of the orbital plane, respectiv ely . If we include the right hand terms of (1), we get the following equations of motion: ¨ r = r ˙ θ 2 − µ E r 2 + ( ~ a pert urb ) r (5a) ¨ θ = 1 r  − 2 ˙ r ˙ θ + ( ~ a pert urb ) θ  (5b) W e no w introduce three assumptions that regard the pertur- bation acceleration terms. First, we ignore all perturbations except for atmospheric drag (i.e., ~ a pert urb ≈ ~ a at mdr a g ). While differential drag takes advantage of the fact that atmospheric drag can secularly affect the in-plane size and shape of an orbit (i.e., its semi-major axis and eccentricity), the other dominant perturbations mainly cause secular changes that are out-of-plane, which we do not ha ve control over . Thus, we omit those perturbations in our simplified simulation model. Second, we assume near-circular orbits where the mag- nitude of the satellite’ s v elocity vector in the tangential direction of the orbital plane is significantly larger than in the radial direction. Since atmospheric drag is antiparallel to the the velocity vector , we ignore the radial acceleration component of the drag perturbation. That is, ( ~ a at mdr a g ) r ≈ 0. W e are left with the following approximated equations of motion: ¨ r = r ˙ θ 2 − µ E r 2 (6a) ¨ θ = 1 r  − 2 ˙ r ˙ θ + ( ~ a at mdr a g ) θ  (6b) Finally , we must approximate the Earth’ s angular velocity in the coordinate frame of the orbital plane so that we may estimate the satellite’ s velocity with respect to the atmosphere (3). Based on the inclination ( φ | 0 ≤ φ ≤ 180 ◦ ) of a satellite’ s orbit, we only consider the component of the Earth’ s angular velocity that is about the normal axis of the orbital plane. W e assume that both the angular velocity of the Earth about its axis and the inclination of the orbit are constant. ~ v rel = r ˙ θ e θ −  ¯ ω E cos ( ¯ φ ) e N × r e r  (7a) = r  ˙ θ − ¯ ω E cos ( ¯ φ )  e θ (7b) Thus, for a near -polar orbit ( φ ≈ 90 ◦ ), the relati ve speed of the satellite is essentially the tangential speed of the satellite (i.e., the speed of the atmosphere is negligible). For a prograde, equatorial orbit ( φ = 0), the relative speed of the satellite is slo wer since the satellite’ s motion is parallel with the velocity of the atmosphere. For a retrograde, equatorial orbit ( φ = 180 ◦ ), the relative speed is faster since the satellite’ s motion is antiparallel with the atmospheric velocity vector . T o summarize, we use the following equations to simulate the approximated orbital motion of a satellite ( ω : = ˙ θ ): ¨ r = r ω 2 − µ E r 2 (8a) ¨ θ = 1 r  − 2 ˙ r ω − 1 2 ¯ C D ¯ m ρ ( r ) | ~ v rel ( r , ω ) | 2 A  (8b) where the notation ρ ( r ) and ~ v rel ( r , ω ) is used to denote that the atmospheric density and relati ve velocity terms are dependent on the satellite’ s radius and angular velocity . D. Appr oximate Discr ete-T ime Model W e introduce a discrete-time dynamical model of the i t h satellite for use in our optimization problem: r i ( k + 1 ) = r i ( k ) + ∆ t · S R ( r i ( k ) , ω i ( k )) · u i ( k ) (9a) ω i ( k + 1 ) = ω i ( k ) + ∆ t · S Ω ( r i ( k ) , ω i ( k )) · u i ( k ) (9b) θ i ( k + 1 ) = θ i ( k ) + ∆ t · ω i ( k ) (9c) + 1 2 ∆ t 2 · S Ω ( r i ( k ) , ω i ( k )) · u i ( k ) Note that our control input is the cross-sectional surface area of the satellite (i.e., u : = A ). The values for S R ( · ) and S Ω ( · ) describe how the impact of the input changes depending on the current state of the satel- lite. For example, for any given cross-sectional surface area, a satellite will experience greater atmospheric drag force at lower altitudes than at higher altitudes. This relationship is captured using the Gaussian variation of parameters (V OP) form of the equations of motion. These equations are used to approximate the rates of change of the time-varying elements in the solution for the unperturbed, two-body system due to small perturbing forces. V allado [15] shows that the a verage rate of change in the semi-major axis of an orbit and the angular speed of the satellite can be expressed in terms of the atmospheric drag perturbation. By applying V allado’ s results to a near-circular orbit, we find the following approximate relationships: S R ( r , ω ) = − ¯ C D ¯ m ρ ( r ) | ~ v rel ( r , ω ) | 2 s r 3 µ E (10a) S Ω ( r , ω ) = 3 2 ¯ C D ¯ m ρ ( r ) | ~ v rel ( r , ω ) | 2 1 r (10b) I I I . O P T I M I Z A T I O N P RO B L E M Our goal is to spread out an initial cluster of satellites in low Earth orbit so that there is equal spacing between each satellite of the shared orbital plane. W e would like to complete this constellation formation maneuver in a fixed number of days while maximizing the operational lifetime of the constellation. The operational lifetime can be defined as the total number of days that all of the satellites remain in orbit. Since atmospheric density increases exponentially as the altitude decreases, a satellite under atmospheric drag experiences very rapid orbital decay as its altitude drops. Thus, our objecti ve is to minimize the drop in altitude of the constellation, which we achie ve by maximizing the altitude of the lowest satellite in the constellation at the final time step T of the optimization problem: maximize U min i = 1 ,..., N r i ( T ) (11) Note that if we ignore oblateness and assume a spherical Earth, maximizing the altitude of a satellite is equiv alent to maximizing the magnitude of its radius. The decision variables are contained in the vector U = [ u 1 ( 0 ) , . . . , u 1 ( T − 1 ) , . . . , u N ( 0 ) , . . . , u N ( T − 1 )] T . So, for a constellation of N satellites and a total of T time steps in which the problem is feasible, there are N T decision variables. T o achieve equal angular spacing of the satellites at the desired final time step T, we use the following inequality constraint in our optimization problem: k D · θ ( T ) − ∆ d es k ∞ ≤ ε θ (12) where θ ( T ) =  θ 1 ( T ) , θ 2 ( T ) , θ 3 ( T ) , . . . , θ N ( T )  T is the angular position state vector at time T . As a conv ention used throughout this paper, we define θ 1 ( t ) as the angular distance trav eled by an arbitrarily designated “lead” satellite in the orbital plane. W e use the same vector notation for angular velocity ω ( T ) and radius r ( T ) . The matrix D is defined as D : =        1 -1 1 -1 . . . . . . 1 -1 -1 1        ∈ R N × N (13) so that D · θ ( T ) =        θ 1 ( T ) − θ 2 ( T ) θ 2 ( T ) − θ 3 ( T ) . . . θ N − 1 ( T ) − θ N ( T ) θ N ( T ) − θ 1 ( T )        ∈ R N × 1 (14) represents the angular separation between adjacent pairs of satellites at the final time step T. W e define ∆ d es : =  2 π N , 2 π N , . . . , 2 π N , - 2 π N ( N − 1 )  T ∈ R N × 1 as the v ector containing the desired angular spacings between each adjacent pair of satellites so that the entire constellation is equally spaced. The dif ference D · θ ( T ) − ∆ d es results in a vector containing the angular spacing errors between adjacent pairs. By constraining the maximum of the absolute values of these errors at time T to be less than or equal to some angular position error tolerance ε θ , we may achieve approximately equal spacing. Similarly , by constraining the angular velocities of adja- cent satellites to be effecti vely zero, we can ensure that the constellation will tend to remain equally spaced in the future and not just at that instance: k D · ω ( T ) k ∞ ≤ ε ω (15) where ε ω is an angular velocity error tolerance close to zero. W e also impose input constraints since our control au- thority is limited by the actual physical dimensions of the satellite: U min ≤ U ≤ U max (16) Giv en the initial state vectors r ( 0 ) , ω ( 0 ) , and θ ( 0 ) , we now summarize the optimization problem here: maximize U min i = 1 ,..., N r i ( T ) subject to k D · θ ( T ) − ∆ d es k ∞ ≤ ε θ k D · ω ( T ) k ∞ ≤ ε ω U min ≤ U ≤ U max (17) which can be restated in the following form by introducing an extra decision variable t : minimize U , t t subject to - r ( T ) ≤ t · 1 N × 1 k D · θ ( T ) − ∆ d es k ∞ ≤ ε θ k D · ω ( T ) k ∞ ≤ ε ω U min ≤ U ≤ U max . (18) W e note that r ( T ) , ω ( T ) , and θ ( T ) do not depend linearly on the input U . T o obtain a linear program, we first precompute reference trajectories ¯ r i ( · ) and ¯ ω i ( · ) by using equations (9) and (10) with the conserv ativ e assumption that each satellite is under minimum drag input until final time step T (i.e., U = U min ). W e then use the following relationships: r i ( k + 1 ) = r i ( k ) + ∆ t · S R ( ¯ r i ( k ) , ¯ ω i ( k )) · u i ( k ) (19a) ω i ( k + 1 ) = ω i ( k ) + ∆ t · S Ω ( ¯ r i ( k ) , ¯ ω i ( k )) · u i ( k ) (19b) θ i ( k + 1 ) = θ i ( k ) + ∆ t · ω i ( k ) (19c) + 1 2 ∆ t 2 · S Ω ( ¯ r i ( k ) , ¯ ω i ( k )) · u i ( k ) where we hav e substituted the reference trajectories in S R ( · ) and S Ω ( · ) so that the equations are linear but time-varying. W e estimate r ( T ) , ω ( T ) , and θ ( T ) from: r i ( T ) = r i ( 0 ) + ∆ t · T − 1 ∑ k = 0  S R ( ¯ r i ( k ) , ¯ ω i ( k )) · u i ( k )  (20a) ω i ( T ) = ω i ( 0 ) + ∆ t · T − 1 ∑ k = 0 n S Ω ( ¯ r i ( k ) , ¯ ω i ( k )) · u i ( k ) o (20b) θ i ( T ) = θ i ( 0 ) + ∆ t · T − 1 ∑ k = 0 ω i ( k ) + 1 2 ∆ t 2 · T − 1 ∑ k = 0 n S Ω ( ¯ r i ( k ) , ¯ ω i ( k )) · u i ( k ) o = θ i ( 0 ) + ∆ t T ω i ( 0 ) (20c) + ∆ t 2 · T − 1 ∑ k = 0 ( T − k − 1 2 ) n S Ω ( ¯ r i ( k ) , ¯ ω i ( k )) · u i ( k ) o which can be expressed in matrix form: r ( T ) = r ( 0 ) + ∆ t · ¯ S R · U (21a) ω ( T ) = ω ( 0 ) + ∆ t · ¯ S Ω · U (21b) θ ( T ) = θ ( 0 ) + ∆ t · T · ω ( 0 ) + ∆ t 2 · ¯ S α · U (21c) where r ( 0 ) , ω ( 0 ) and θ ( 0 ) ∈ R N × 1 are the initial state vectors and ¯ S R , ¯ S Ω are large matrices of the form: ¯ S R =    ¯ S R 1 . . . ¯ S R N    , ¯ S Ω =    ¯ S Ω 1 . . . ¯ S Ω N    ∈ R N × ( N · T ) (22) consisting of the following row vectors along the diagonal for i = 1 , . . . , N : ¯ S R i = [ S R ( ¯ r i ( 0 ) , ¯ ω i ( 0 )) , . . . , S R ( ¯ r i ( T − 1 ) , ¯ ω i ( T − 1 ))] (23a) ¯ S Ω i = [ S Ω ( ¯ r i ( 0 ) , ¯ ω i ( 0 )) , . . . , S Ω ( ¯ r i ( T − 1 ) , ¯ ω i ( T − 1 ))] (23b) The ¯ S α matrix is of the same form where the diagonal ro w vectors are found by calculating: ¯ S α i =  ( T − 1 2 ) · 1 1 × T − [ 0 , . . . , ( T − 1 )]  ◦ ¯ S Ω i (23c) where ( ◦ ) is the element-wise product of the two ro w vectors. The ¯ S R , ¯ S Ω , and ¯ S α matrices are precomputed prior to solving the program, which we can now express in standard form: minimize x f T x subject to Ax ≤ b (24) where x = [ U , t ] T and f T = [ 0 1 × ( N · T ) , 1 ] . The matrices used to form the inequality constraints are: A =           - ∆ t · ¯ S R , - 1 N × 1 ∆ t 2 · D · ¯ S α , 0 N × 1 - ∆ t 2 · D · ¯ S α , 0 N × 1 ∆ t · D · ¯ S Ω , 0 N × 1 - ∆ t · D · ¯ S Ω , 0 N × 1 I ( N · T ) × ( N · T ) , 0 ( N · T ) × 1 - I ( N · T ) × ( N · T ) , 0 ( N · T ) × 1           (25) and b =           r ( 0 ) ε θ · 1 N × 1 − D · [ θ ( 0 ) + ∆ t · T · ω ( 0 )] + ∆ d es ε θ · 1 N × 1 + D · [ θ ( 0 ) + ∆ t · T · ω ( 0 )] − ∆ d es ε ω · 1 N × 1 − D · ω ( 0 ) ε ω · 1 N × 1 + D · ω ( 0 ) u max · 1 ( N · T ) × 1 - u min · 1 ( N · T ) × 1           (26) I V . S I M U L A T I O N R E S U LT S The simulation results in this section are based on N = 105 satellites beginning at the same initial states, corresponding to a Sun-synchronous, circular orbit at an altitude of 475 km: θ i ( 0 ) = 0, r i ( 0 ) = r 0 and ω i ( 0 ) = q µ E r i ( 0 ) 3 for all i = 1 , . . . , N . In our linear program, we set the angular separation and velocity error tolerances at ε θ = 0 . 1 degrees and ε ω = 1 e -18 rad/s, respectively . The solution to the program corresponds to optimal drag area commands (i.e., desired cross-sectional surface areas) that are sent simultaneously to all satellites once every 24 hours. The drag area commands are allowed to take on any v alue within a continuous range between the minimum and maximum possible surface areas. A. Open-loop versus F eedback Starting at the initial conditions, we determine that the linear program is feasible with a horizon of T = 71 days. W e apply the optimal input commands in open-loop and find that after 71 days, the angular spacings between adjacent satellites tend to conv erge tow ards the desired value of 2 π N (see Fig. 1). Howe ver , due to the approximations made in the linear model (19), none of the angular spacings satisfy the angular spacing error tolerance that we defined. Furthermore, although the program expects a maximum altitude drop of 10.79 km at the end of 71 days, the open-loop simulation results in a worse 11.64 km altitude drop. Fig. 2 shows how the solution to our initial optimization problem is to utilize the whole range of input values, from minimum to maximum, and to giv e each satellite a different input v alue at time step k = 0 based on its arbitrary number- ing (i.e., the 1 st satellite receiv es a maximum drag command while the N t h satellite receives a minimum drag command). By the final time step T , the input value for any particular satellite is “flipped” in intensity compared to its initial value. For example, the input for the 1 st satellite changes from an initial maximum drag command to a final minimum drag command. All the while, the “middle” satellite receives a relativ ely steady input value at all steps. The result is that the satellites which are arbitrarily assigned lo w or high numbered positions (i.e., 1, 2, 3, ... or ... ( N − 2 ) , ( N − 1 ) , N ) e xperience larger changes in input than the satellites placed to wards the middle of the constellation position numbering scheme. Due to the dissatisfying open-loop performance, we then use a model predictiv e control (MPC) approach where the linear program is solved at the beginning of each time step Fig. 1: Optimal input commands applied in open-loop result in final angular spacing values that fail to land within the designed error tolerance thresholds (represented by dashed lines). Note that the angular spacing between the N t h and 1 st satellites ( θ N − θ 1 ), which should con verge to - 2 π N ( N − 1 ) , is not included in this figure. Fig. 2: Optimal input commands for all N satellites applied in open-loop each day until horizon T = 71 days. but only the first set of control inputs in the sequence is applied to the satellites. At the first time step, we solve for the horizon T . At each subsequent time step, the problem is reformulated and solv ed again but with a shrinking horizon (i.e., at k = 0, horizon = T , at k = 1, horizon = T − 1). By solving the program at each time step with updated states, we are able to correct for prediction errors and the effect of un-modeled perturbations. As can be seen in Fig. 3, when the optimal inputs are applied with feedback, an equally-spaced constellation is formed within 71 days where all the angular spacings satisfy the designed angular separation error tolerance. Also, compared to the 11.64 km altitude loss in the open-loop simulation, Fig. 4 shows that all the satellites under feedback control con verge in altitude and angular velocity at the expense of losing only 10.71 km in altitude (i.e., we conserv e 1 km of altitude), which is very close to the 10.79 km predicted by the program at time step k = 0. This 10.71 km drop in altitude compares to a 2.84 km drop under constant minimum drag and a 19.88 km drop under maximum drag, for the same number of days. Fig. 4 also shows how the altitude of each satellite is varied over time, resulting in an inv ersely proportional change in angular velocity . This difference in angular velocity between pairs of satellites allows the controller to adjust the angular spacings to the desired values. Fig. 5 shows the input area commands that are computed and applied at the beginning of each day . Fig. 3: Angular spacing between adjacent pairs of satellites change from 0 ◦ to 360 ◦ N ± ε θ . Note that the angular spacing between the N t h and 1 st satellites ( θ N − θ 1 ), which should con verge to - 2 π N ( N − 1 ) , is also not included in this figure. B. T rading Constellation Acquisition T ime for Lifetime W e increase the horizon length T (i.e., number of days for the satellites to con ver ge to the desired constellation) to de- termine the effect on the overall lifetime of the constellation. Fig. 6 sho ws that as we increase the horizon length from 71 days to 98 days, we can reduce the altitude drop by 2.43 km, effecti vely extending the lifetime of the constellation. Howe ver , increasing the horizon beyond 98 days does not result in further improv ement. C. Maintaining Constellation thr oughout Lifetime Once the equally-spaced constellation is achie ved, the initial acquisition phase is complete and the satellites enter the operational mode where they are allowed to “drift” in a minimum-drag attitude configuration. Howev er, when the spacing between any adjacent pairs reaches abov e a certain threshold, we apply the optimal control strategy again, albeit with a much shorter horizon. In Fig. 7 we sho w that our approach is successful in maintaining the angular separations throughout the operational phase. W e also observe that the Fig. 4: At the horizon T , the altitudes (and angular velocities) of the satellites con verge, ensuring that the constellation will remain in the desired configuration with minimal control effort for the remainder of the constellation lifetime. Fig. 5: The average lev el of actuation is high both in the beginning and tow ards the end of optimal constellation separation phase. orbital motion of the constellation is relati vely “smooth” in the operational phase compared to the initial acquisition phase (see Fig. 8). Furthermore, in Fig. 9 we see that the lev el of actuation required to maintain the constellation is significantly less than that required in the acquisition phase. W e arbitrarily define the constellation operational lifetime to be the total number of days that none of the satellites drops to an altitude of 200 km or less, where spacecraft orbits decay rapidly . The total constellation lifetime is 1,059 days (2.90 years) under this optimal angular separation control strategy . As a comparison, a satellite under constant minimum or maximum drag would ha ve a lifetime of approximately 1,410 days (3.86 years) or 232 days, respectiv ely . Fig. 6: Smallest altitude drop of 8.28 km is achie ved when horizon length is 98 days. Fig. 7: Angular separation is maintained for the lifetime of constellation. Note that the angular dif ference between the N t h and 1 st satellite reaches a value of - 360 ◦ N ( N − 1 ) ± ε θ . V . C O N C L U S I O N S Although the orbital dynamics are nonlinear , we found that both the altitude and angular velocity of a satellite controlled by differential drag can be approximated as linear over the relativ ely small operating range of a singe day . Thus, the solution from our linear program, e ven when applied in open- loop, pro vides reasonable performance in forming the equally spaced constellation. T o compensate for the prediction error caused by model-process mismatch and improve controller performance, we lev erage the feedback mechanism provided by the shrinking-horizon MPC approach. W ith feedback, we are able to achieve an equally-spaced constellation that satisfies design tolerances and av oids unnecessary control action that reduces operational lifetime. W e also observed that we can increase the operational lifetime of the con- stellation by allowing it to form over a longer time frame. Fig. 8: Once equally spaced, the orbital motion of the constellation is relativ ely “smooth” compared to the initial, optimal constellation separation phase. Fig. 9: Once equally spaced, relatively minimal actuation is required to maintain the angular separation for the duration of the constellation operational lifetime. Howe ver , the tradeof f can only be made until a certain point at which increasing the horizon further results in decreased lifetime. Finally , we show that the constellation can be maintained throughout its lifetime by applying the same optimal control strategy when the angular spacing errors drift abov e a designed threshold value. AC K N OW L E D G M E N T Thanks to Professor Simone D’Amico of the Department of Aeronautics and Astronautics at Stanford University for introducing this problem in the AA279A: Space Mechanics course of Winter 2017. The authors gratefully acknowledge support from the National Science Foundation under grant ECCS-1405413. Andrew Packard ackno wledges the generous support from the F ANUC Corporation. R E F E R E N C E S [1] J. Puig-Suari, C. Turner , and W . Ahlgren, Development of the standard CubeSat deployer and a Cubesat class Picosatellite, IEEE Aerospace Conference Proceedings, vol. 1, 2001. [2] K. W oellert, P . Ehrenfreund, A.J. Ricco, H. Hertzfeld, Cubesats: cost-effecti ve science and technology platforms for emerging and dev eloping nations, Adv ances in Space Research, vol. 47, pp. 663?684, 2011. [3] S. Bandyopadhyay , R. Foust, G.P . Subramanian, S.J. Chung, and F .Y . Hadaegh, Review of Formation Flying and Constellation Missions Using Nanosatellites, Journal of Spacecraft and Rockets, pp. 1-12, 2016. [4] D. Selva, D. Krejci, A survey and assessment of the capabilities of Cubesats for Earth observation, Acta Astronautica, vol. 74, pp. 50-68, 2012. [5] Y . Hu and V . Li, Satellite-based Internet: a tutorial, IEEE Communi- cations Magazine, vol. 39, no. 3, pp. 154-162, 2001. [6] F . Khan, Mobile Internet from the Heavens, http://arxiv .org/abs/1508.02383. [7] C.L. Leonard, W .M. Hollister, and E.V . Bergmann, Orbital Forma- tionkeeping with Differential Drag, Journal of Guidance, Control, and Dynamics, vol. 12 no. 1, pp.108-113, 1989. [8] B. S. Kumar , A. Ng, K. Y oshihara, and A. De Ruiter, Differential Drag as a Means of Spacecraft Formation Control, IEEE Transactions on Aerospace and Electronic Systems, vol. 47, no. 2, pp. 1125-1135, 2011. [9] B.D. Mills, Satellite Paradox, American Journal of Physics, vol. 2, no. 2, pp. 115-117, 1959. [10] L. Blitzer , Satellite Orbit Paradox: A General V iew , American Journal of Physics, vol. 39, pp. 882-886, 1971. [11] A. Li and J. Mason, Optimal Utility of Satellite Constellation Sep- aration with Differential Drag, AIAA/AAS Astrodynamics Specialist Conference, 2014. [12] C. Foster, H. Hallam, and J. Mason, Orbit Determination and Differential-drag Control of Planet Labs Cubesat Constellations, AIAA/AAS Astrodynamics Specialist Conference, 2015. [13] R.R. Bate, D.D. Mueller, J.E. White, Fundamentals of Astrodynamics. New Y ork: Dover , 1971, ch. 1. [14] O. Montenbruck and E. Gill, Satellite Orbits. New Y ork: Springer, 2005, ch. 3. [15] D.A. V allado, Fundamentals of Astrodynamics and Applications. Cal- ifornia: Microcosm, 2013, ch. 9.