A network flow approach for constellation planning

A graph-based network flow approach to constellation planning is presented that respects spacecraft resources and mission constraints to arrive at a preliminary operational schedule. The method discretizes task fulfillment opportunity windows into graph nodes and then adds edges based on feasibility of transition between these opportunities, either by slew maneuvers or extended-time task fulfillment during downlink operations. Intelligent graph pruning is applied to improve computation performance and provide a versatile planning capability with a minimal impact on the resulting solution. A network flow formulation allows for a simultaneous search for spacecraft plans using a Mixed Integer Linear Program (MILP). The framework enables the definition of appropriate mission constraints to ensure viability and mission efficiency among the constellation of satellites. The planning technique is demonstrated for a mission of 100 satellites to illustrate its capability and performance within this proliferated regime.


INTRODUCTION
Satellite constellations have been conducting operations in space since the early 1960s when the CORONA mission became operational [1] .Since then, a number of additional constellations have joined this exclusive club.These include new Earth observing constellations such as Maxar' s WorldView [2] , NASA' s A-train [3] , and Planet Labs' Doves [4] , to name just a few.While many of the individual satellites are shrinking in mass, volume, and cost [5] , the constellations continue to grow in numbers, capability, and sophistication, requiring increased levels of automation [6] to properly orchestrate in accomplishing the mission objectives [4,7] .As new mega-constellations are designed and deployed, several new challenges are realized [8,9] .Some of these challenges deal specifically with general geometric configuration and orbital design [10,11] , with some research recommending methods to address the challenge [12] .Additionally, and more importantly to the content of this paper, is that as the number of satellites within a constellation increases, the operational planning problem complexity compounds due to the spatial and temporal dependencies present within the constellation' s operational environment, and the limited resources onboard each vehicle [13] .
The Earth imaging geospatial intelligence (GEOINT) mission, addressed in this paper, falls into the broader category of an assignment problem where the goal is to optimally assign tasks to agents [14] .It can be further classified into a particular group using a well-defined taxonomy developed by Gerkey and Mataric that describes different types of tasks and the agents capable of fulfilling those tasks [15] .The following three axes of differentiation are used for that classification: 1.The number of tasks a robot/agent can perform -a single task (ST) or multiple-tasks (MT) at a time.2. The number of robots/agents required by a task in order to fulfill it -single robot (SR) tasks, or multi-robot (MR) tasks.3. The time of allocation -instantaneous allocation (IA) or scheduling robots/agents over an extended period of time (TA).
Korsah et al. extended this concept by accounting for dependencies that often exist between the tasks and the robots/agents executing them [16] .These dependencies are the following: 1.No Dependency (ND) -no dependency exists between tasks or agents.2. In-schedule Dependency (ID) -tasks fulfillment decisions impact other tasks fulfillment opportunities.
3. Cross-schedule Dependency (XD) -agents must coordinate schedules and collaborate directly for task fulfillment.4. Complex Dependency (CD) -task fulfillment utility depends on the schedules of other agents in the system in a manner that is determined by the particular task decomposition chosen.
A note of particular importance, made by Gerkey and Mataric, is that any problem beyond the [ST-SR-IA] is strongly NP-hard and thus precludes enumerative solutions due to computational complexity and the associated run time [15] .In this paper, the assignment problem is that of tasking multiple satellite agents to optimally fulfill the GEOINT mission, which requires efficient image collection and data delivery to the ground.Using the taxonomies defined, this problem falls into the ID category due to having dependencies on the other tasks a particular agent is performing.Furthermore, each vehicle is only capable of performing a single task at a time, and that task only requires a single satellite to fulfill it, thus placing this problem into the [ST-SR-TA] category.The full categorization is then ID[ST-SR-TA] that, while not being the most complex problem, is still known to be NP-hard [16] .
The GEOINT design reference mission (DRM) constellation scheduling problem is of special interest at this time to the space community and deserves additional research attention due to the complexity of formulating and efficiently solving it for mission-realistic scenarios.The primary objective of this mission type is to generate an efficient, operationally-viable schedule that fits within all specified constraints of the mission.Several approaches have been presented to address a portion of this DRM.A host of work has focused on arriving at feasible, and potentially optimal, solutions to the image collection problem [17][18][19][20][21][22][23] .Similarly, the data relay portion of the problem has also been investigated by several researchers and they provide robust solutions to scheduling the transfer of data [24][25][26][27][28][29][30] .However, in both the image collecting maximization and data relay scheduling problems, the papers mentioned address them independently.In reality, the image collection and downlink problems are directly coupled and should be solved simultaneously to yield a reliable and accurate solution for the constellation.
Some works have considered both the image collection and data relay problems within the same framework and provided various solutions to solve them [31][32][33] .However, Hu et al. dedicated most of their paper to evaluating a branch and price technique for finding a solution instead of also providing an in-depth explanation of the problem itself [32] .Additionally, Hu et al. simulated a small LEO constellation of only 3 satellites and did not evaluate the solution methodology in the proliferated LEO context that the industry is beginning to witness today [34] .Similarly, Peng et al.only formulated their planning problem for a single satellite and neglected its application to a multi-agent cooperative system [33] .It is the opinion of the authors of this paper that a truly applicable solution to the GEOINT constellation scheduling problem must be capable of simultaneously planning operations for at least 100 Earth imaging satellites and delivering a robust solution for each of the individual spacecraft elements.The paper by Augenstein et al. addressed many of the requirements for a mission-relevant constellation planning technique and was of particular interest to the authors.Augenstein et al. formulated the problem for many cooperative agents and provided a solution to the coupled image collection and data downlink problem.However, they required a separation of the coupled problem and relied on heuristics in order to solve it within the required timeframe for their application.Doing so, enabled a rapid solve time but sacrificed a globally-optimal solution within the discretized space.To address these issues, the authors of this paper leveraged the work from Augenstein et al. and augmented it with a novel approach using network flow theory, which provided an alternative method for specifying mission objectives, vehicle constraints, and additional methods for improving the solve time.This work directly confronts the challenges posed by the constellation planning problem and presents an optimal planning technique for the GEOINT DRM that enables the simultaneous coordination of both image collection and scheduling of ground stations for data downlink, across the cooperative constellation of space vehicles.
The network flow technique, described in this paper, discretizes the possible task fulfillment opportunities for each satellite and represents available transitions between tasks in graphical form.This graphical formulation enables a simultaneous search for each satellite' s path through fulfillment of image collection and downlink tasks.The application of network flow theory to this constellation problem is one of the central contributions made within this paper.Flow constraints allow for coordination across the constellation while resource constraints ensure a mission-realistic solution.The network flow formulation is used to create an operational schedule using a Mixed-Integer Linear Program (MILP) and can be used to solve the constellation planning problem in a mission-relevant way.Instantaneous tasks (e.g., image collection) are represented as nodes in the graph while extended-time duration task fulfillment opportunities (e.g., downlinking of data or slewing the satellite) are modeled as edges in the graph.Fundamental constraints such as lighting and on-board satellite memory are also considered.This blending of task generation, and fulfillment, within an operationallyconstrained environment builds on previously introduced planning methods in the literature.
The remainder of the paper will proceed with a brief background on constellation planning methods and algorithms in Section 2. This is followed by a detailed decomposition of the problem into a graph-based representation in Section 3 and constructive development of a network flow optimization solution in Section 4. Finally, a mission scenario of 100 sun-synchronous imaging satellites is analyzed using the network flow approach to illustrate realistic mission performance in Section 6.The paper ends with concluding remarks in Section 7.

BACKGROUND AND RELATED WORK
As mentioned above, the purpose of a constellation planning system, within the DRM context, is to maximize prioritized image collection and data downlink for mission task fulfillment while respecting each satellite' s capabilities and available resources.The satellite constellation and ground resources must be treated holistically to ensure that satellites do not duplicate efforts in imaging and that only a single satellite uses a particular ground terminal at any given time for downlink of mission data.These problem elements are directly coupled since image data collection must occur prior to that data being sent to the ground via downlink.Additionally, the slew agility and onboard memory capacity of individual satellites must be considered to ensure an executable and mission-relevant plan.Finally, while image collection tasks can be performed nearly instantaneously, the downlink tasks require extended periods of time to fulfill.
Augenstein et al. discretize both downlink and image collection opportunity time windows to initialize the search space considered in arriving at a solution for this coupled problem [31] .While the solutions to the discretized problem are inherently sub-optimal to the granularity of the discretization, the discretization enables the use of a MILP formulation, making the combinatorial problem tractable.Task fulfillment instances allow the planner to evaluate specific decisions within a constrained and finite environment.
To increase applicability to realistic mission scenarios, several constraints are identified by Augenstein et al.
The satellites can only perform a single task at a time, either image collection or downlink, and they are constrained with a finite slew rate, meaning that pointing maneuvers can only be conducted up to a specified rate.Satellites also have a fixed amount of onboard memory that can be consumed by image data storage [31] .
Ground stations can only communicate with a single satellite at a time, requiring a deconflicted schedule for ensuring data delivery to the ground.An additional operational constraint is applied in that every 3 orbit revolutions, each satellite contacts the ground for at least 3 contiguous minutes.This helps maintain up-to-date ground knowledge about the health and safety of each vehicle within the constellation [31] .
The combinatorial optimization problem addressed by Augenstein et al., is defined as a MILP with constraints and objectives applied globally rather than using the original formulation of a DAG for each satellite.This MILP formulation can be solved to arrive at an optimal solution within the discretized space.However, the authors point out that the run-time performance of the formulated MILP is unacceptable to their mission timeline since they require multiple revisions to the plan based on immediate mission needs and potential late tasking.To address this problem they separate the immediate coupling that exists between image collection opportunities and downlink opportunities.Separating the imaging and downlink planning decouples the variables present in the MILP but requires a method to appropriately allocated downlink time based on estimated image collection and current onboard memory state of each vehicle.To accomplish this, they provide a heuristic that estimates the amount of priority-weighted image data likely onboard the satellite as a function of time and the number of priority-weighted image collection opportunities a satellite would encounter during specific time periods.The authors demonstrate that with these two pieces of information it is possible to reasonably allocate downlink time among the constellation of vehicles and then also schedule the additional image collection opportunities.This modification yields dramatically reduced solve times and enables the efficient, yet sub-optimal, scheduling of the TerraBella constellation of satellites [31] .
The approach developed herein provides several contributions to the planning framework established by Augenstein et al.For example, the discretized task opportunities are organized into a directed acyclic graph (DAG) with network flow constraints used to allow for simultaneous search for each satellite' s plan.Utilizing the DAG formulation, transition costs can be considered instead of solely scoring the imaging opportunities.Furthermore, the amount of data to be downlinked can directly be associated with the time spent along edges connecting downlink nodes rather than only scoring the instantaneous node as in Augenstein et al.By doing so, the edge transitions, and their associated costs, more effectively represent the actual operations to be conducted on orbit.This method is thus more realistic and is more intuitive for an operations team to understand.
Additionally, the flow formulation enables the definition of concise application of mission constraints to avoid duplication of task assignment and does not rely on heuristics to solve.The intent of this research is to formu-

Discrete time index
late the problem in such a way that an optimal solution is possible to obtain within a mission-realistic timeline.While the MILP solver does use various heuristics during processing, mission-specific heuristics were not developed.Instead, the formulation discusses a graph pruning concept to simplify the search space without sacrificing optimality.As the DAG can become unwieldy with large numbers of edges, another contribution is the technique of trimming unnecessary edges with minimal effect on the resulting solution.A minor contribution is the explicit definition of image collection constraints such as allowable target collection geometries, target lighting, etc.A complete reference to the notation used in this paper is provided in Table 1.

GRAPH-BASED CONSTELLATION PLANNING
The network flow constellation planning approach taken herein can be broken into several distinct steps, as illustrated in Figure 1.It begins by computing the access windows of satellites to targets and ground station antennas.These access windows are defined by start and end times and can be discretized over that period of time.The specific instances in time become nodes within the graph.Nodes are then connected by edges if the transition is feasible, meaning that the required slew rate is below the satellite' s allowable slew rate.Once the nodes and edges have been established, a utility is associated with each edge based upon the transition cost (e.g., slew rate) and the score (e.g., downlink value) or imaging task score at the node.The various constraints on task  The minimum and maximum azimuth angle from the target to the satellite min/max elevation The minimum and maximum elevation from the target to the satellite min/max sun elevation The minimum and maximum sun elevation angle from the target start/end time Temporal limits on the valid time window in which the satellite can perform the task fulfillment and data collection are then associated with the edges.The final step is to solve for the paths through the graph which optimize the utility, subject to the operational constraints applied.The subsections that follow provide precise task constraint definition and greater detail on each of the steps taken during constellation planning.

Task definition
The satellites, within the constellation, will be commanded to perform two distinct types of tasks.The first is an imaging task and the second a downlink task.The imaging task corresponds to collecting an image of a particular target and the downlink task corresponds to downlinking image data to a ground station.Both the imaging and downlink tasks have temporal and spatial constraints dictating how and when the task can be performed.The imaging task may have an additional lighting constraint.Each of these constraints is summarized in Table 2.The spatial and lighting conditions can be transformed into temporal conditions for each satellite by considering the satellite trajectory, which defines the position of the satellite as a function of time.Each constraint is now discussed by assuming that a time value and satellite position are specified.
The temporal constraint satisfaction is determined by evaluating whether or not the opportunity time occurs within the allowable execution window.Temporal constraints allow the system to set limits on when a target is available for collection due to seasonal or weather concerns, or when a ground system outage occurs for a downlink antenna.
The spatial constraints depend upon the position of the satellite.These constraints are defined in terms of elevation and azimuth angles and use frames centered at the imaging target or ground antenna.Generally, the elevation angle constraints for ground antennas are set to ensure line-of-sight access to the satellite from the antenna, whereas the elevation constraints for an imaging target ensure high overhead imaging to ensure quality data.The azimuth angle constraint controls which side is available to the satellite and can be used to avoid target or antenna occlusion from neighboring buildings or nearby mountainous terrain.
The calculations for the elevation and azimuth angles are performed in the East-North-Zenith (ENZ) frame as depicted in Figure 2, and denoted as  in this paper.The relative position of the satellite with respect to the target in this frame is given as ( The elevation angle of the satellite relative to the target is calculated as ( The azimuth angle, of the satellite relative to the target, is calculated in a similar manner as described by (3) In addition to the temporal and spatial constraints, the imaging task may also include lighting constraints.For example, the simulation in Section 6 constrains all targets to be valid only during the daytime (i.e., sun elevation angle between 0°and 90°).The calculation of the sun elevation angle is done exactly as before by replacing the satellite location with the apparent location of the sun, as in The time in question is considered valid for a particular satellite if the time is within the temporal bounds and the calculated elevation and azimuth angles from the target to the satellite, as well as from the target to the sun, are within the angular bounds at the given time.Given the valid access times, access windows for each task are found using root finding techniques to determine the start and end time of each window.

Graph node generation
The access window of each task forms an essential element to the creation of graph nodes.There are four types of nodes added to the graph.A starting and ending node for each satellite is added to the graph corresponding to the satellite position at the start and end time of the planning horizon.Backbone nodes are added (discussed in Section 5 as means to maintain connectivity of the graph when pruning edges).The fourth node type corresponds to a task (image collect start time or ground station opportunity start time discretized as shown in right-hand side of Figure 2).The nodes are denoted as   ,  = 1, ...,   where   is the total number of nodes in the graph.Each node corresponds to a specific satellite at a specific point in time, either performing a task or in standby.
The task nodes are created from a discretization of the access window for each task with each point of discretization becoming a node in the graph.A more accurate graph can be created with a finer sampling of the valid time windows, but at the cost of increasing   , which adversely affects the graph complexity.Figure 2 illustrates a valid access window for a particular target and satellite flight path.The time window discretization is overlaid on the valid access window to arrive at the initial set of nodes for that particular target and satellite.
An orientation vector,   where  is the node index, is associated with each node except the ending node.For the starting node,   corresponds to the pointing vector of the camera at the starting time.For imaging and downlink nodes,   corresponds to the vector pointing from the satellite to the target or ground antenna at that point in time while the backbone nodes assume a nadir orientation, and the final graph node is not assigned a value for   .Not associating an orientation with the final node effectively allows any node to be the final imaging task assigned for the satellite.
Also associated with node  is a score,   ≥ 0. As each node is associated with a particular satellite at a particular point in time, the score for a particular task can vary based upon the satellite performing the task and the time at which the task is being performed.Establishing viable transitions between these nodes is discussed in the next section.

Graph edge creation
With the nodes established, it is now possible to evaluate the feasibility of the satellite to transition between them via an edge.Associated with each potential edge is a slew angle that corresponds to the change in orientation between nodes.The required slew rate is determined by calculating the required angular change divided bye allowable time between nodes.An edge is only added to the graph if the slew rate is below the allowable slew rate, defined by the agility of the particular satellite.
Given two pointing vectors,   and   corresponding to nodes   and   , the angle between them,    , can be calculated using the relation between the cosine and dot product: Using   and   to denote the times that correspond to nodes   and   , the nominal slew rate is calculated as A transition between nodes is deemed feasible if the nominal slew rate is below the maximum allowed slew rate for a particular satellite.If feasible, the edge is added to the graph.As no orientation angle is associated with the final node, an edge is allowed between all nodes and the final node.This effectively allows any node to transition to the final task guaranteeing at least one valid solution through the graph.

Graph scoring and costs
When performing a graph search, the optimization engine attempts to traverse the graph in a way that maximizes utility for the mission.Utility is defined as the score of performing a task minus the cost, or penalty, required to do so.The globally-optimal plan considers the mission holistically and selects a path through the graph for each satellite that accomplishes this.Each satellite is not necessarily optimal but the complete mission plan is.
For imaging tasks, the score is the parameter used for assessing priority of the collection and acts as a reward for arriving at the end node of the edge (i.e., the worth of the task performed at the end node).The notation   represent the node index at the end of each edge , denoted   .For example, if   connects nodes  and  , then   =  .The notation   is used to represent the score of an edge.Imaging tasks scores are contained on the edge entering the imaging node.
For downlink tasks, the score for task  is calculated as the duration of the edge (  ) multiplied by the data rate () and the arriving edge score (   ), divided by the on-board memory size of the satellite (), or ( For edges that end at a downlink task but do not come from a downlink task for the same ground station the score is 0,   = 0.This is because the graph must enter a downlink start node before realizing the downlink score when traversing the subsequent edges. The cost, or penalty, associated with each edge of the graph is the weighted slew rate necessary to perform the maneuver between neighboring task opportunities (e.g., from node ℎ to node ), with a higher required slew rate resulting in a higher penalty.The utility at node , when transitioning from node ℎ, is given as where the slew cost weighting factor is represented as .The utilities are combined into a single vector as , where   is the total number of edges.

The graph representation
The graph  is now formally defined and contains set of nodes (or vertices), , and edges,  ⊂  × , i.e.  = {,  }.Given two nodes in the graph,  1 ,  2 ∈ , an edge from  1 to  2 implies that the pair ( 1 ,  2 ) is in the edge set, , i.e. ( 1 ,  2 ) ∈ .An example graph is shown in Figure 3 with four nodes and five edges.Recall that each node corresponds to an instant in time and applies to a particular satellite.Thus, given a constellation of   satellites,  consists of   distinct subgraphs.
The planning will make use of the graph' s incidence matrix, , as defined in [35] .Note that the definition of incidence matrix often includes a positive one in both nonzero rows of each column.The incidence matrix as defined herein is also referred to as the coefficient matrix [36] .Given | | =   and | | =   , where | • | denotes the cardinality of a set, the incidence matrix is of dimension   ×   .Each column of  is used to represent an edge.Expressing the  ℎ row of the  ℎ column as    , the incidence matrix can be expressed element-wise as follows: where Edge  originates at node  −1 Edge  ends at node  0 Otherwise (7)   The signed elements of the incidence matrix allow for quick evaluation of the number of incoming and outgoing edges for each node.This will be useful later to generate continuous paths through the graph.One additional matrix that will be used for evaluating the density of the graph is the directed adjacency matrix.The adjacency matrix, denoted , is an   ×   matrix where row  column  will equal 1 if an edge begins at node   and ends at node   , with a value of 0 otherwise.

Simple Example of the Incidence and Adjacency Matrices:
For example, the incidence and adjacency matrices for the graph in Figure 3 are written as Column  of  corresponds to edge   in Figure 3 with the positive element corresponding to the originating node.Row  of the adjacency matrix shows the nodes to which node  can directly connect.With nodes and edges connected and properly quantified for utility, the flow through the graph can now be investigated.

A NETWORK FLOW APPROACH FOR CONSTELLATION PLANNING
Our attention now turns to formulating a network-flow based optimization problem for constellation planning.
The constellation planning problem that includes image acquisition and downlink scheduling will be solved by finding a set of paths within the DAG with the maximal utility (score-penalty), subject to the mission constraints of non-overlapping task assignments and the limitations imposed by each satellite.These include on-board memory and slew agility.Iterative, dynamic programming-based techniques such as Dijkstra and its many derivatives, e.g., [37] , could be used to find a solution to the unconstrained graph search problem for a single satellite.However, a batched solution using a MILP formulation will be used as it can readily incorporate the constraints required for imaging and be used to simultaneous plan for all satellites.This section constructively develops the MILP problem.First, a network flow approach is presented for graph search.A group constraint is then added to ensure that only one of the many nodes corresponding to an imaging task will form part of the solution.A satellite memory constraint is then added to ensure proper downlinking of information to respect onboard data storage capacities.

A Network flow approach to graph search:
A mixed integer formulation of the graph search is developed by using a binary decision variable to represent each edge in the graph, denoted as   ∈ {0 1} for  = 1, ...,   (recall that   = | |).A value   = 1 indicates that edge  is part of the path and a value of 0 indicates that edge  is not.To select one path over another, the utility   discussed in Section 3.4 is associated with each edge.The total utility of the selected paths can be written as the summation of the individual utilities,  1  1 +  2  2 + ... +       =   .
To ensure that the choice of   forms continuous paths from each satellite' s starting node to each satellite' s ending node, a network flow conservation constraint is employed.This constraint requires that the start node have a single exit path, the end node has only one entry path, and all other nodes, between the start and end node, have both a single entry and single exit to conserve flow.Imagine a simplistic network of pipes carrying water.Each node represents a connection point of various pipes and each edge represents the pipes.Special nodes called "sources" can provide water while others called "sinks" can store or consume the water.All other nodes are intermediary and must simply pass out whatever water comes into the node.Returning to the path planning problem, the "source" nodes are where the satellites start and the "sink" nodes are the terminal nodes for each satellite.All intermediary nodes are decision points.See [36] Chapter 10 for a thorough review of network flow problems.
For the graph search problem, the network flow conservation constraint consists of the source providing a unit of flow, the sink accepting a unit of flow, and all other nodes having a balance of flow (i.e., if an incoming edge to the node is selected, exactly one outgoing edge must also be selected).Recall that row  of the incidence matrix will have a "1" in columns corresponding to edges that start at node  and a "−1" where edges terminate at node .Writing   as the  ℎ row of the incidence matrix,    =   implies that there are   more edges originating at node  than terminating at node  (note,   can be negative).Thus, a source node must have   = 1, a sink node   = −1, and an intermediate node   = 0.
Assuming that no two satellites are collocated, each satellite will have a distinct access schedule to each imaging task.The aggregate graph nodes and edges could each be separated into   disjoint subsets, one for each satellite.
Each disjoint subset has its own source and sink node.  is used to denote the incidence matrix corresponding to the access schedule of satellite .Assuming that the source node for each satellite corresponds to the first row of   , the sink node corresponds to the final row, and an abuse of notation is used to write the edge variables corresponding to satellite  as   , the network flow constraint for satellite  can be represented as The aggregate network flow could be written by combining the disjoint components of the incidence matrix as follows.
The optimization problem can then be written to include this constraint as: Simple Example of the Network Flow Constraint: Consider the simple example from Figure 3 with corresponding incidence and adjacency matrices in (8).The network flow constraint in (9) would be reduced to A line-by-line overview of the constraint can help to understand the flow constraint.Line 1 indicates that  1 must be used.This is also obvious from Figure 3 as it is the only edge coming out of the source.Line five states that either  6 be active or  7 , but not both, which is also readily apparent from the figure as these edges lead into the sink.Rows two, three, and four maintain the balance of incoming edges to outgoing edges and ensures the continuous flow through those nodes.

The group constraint
As each satellite has its own disjoint subgraph, each satellite could plan independently.However, multiple nodes in the graph exist for each task.A group constraint is added to provide the requisite coordination, ensuring that only one node for each task will be selected.Each group of nodes is collected into a set where V  ⊂  corresponds to group  and there are a total of   groups (i.e.,  = 1, ...,   ).The set of edges leading into group  is denoted as E  , where Note that this edge subset may include multiple edges corresponding to a single satellite' s disjoint edge set or edges from multiple satellite edge sets.The index set, I  , contains the edge indices corresponding to E  .The constraint for group  can be written in summation form as Combined with the fact that   ∈ {0, 1}, (10) ensures that a maximum of one edge will enter the group.
I  can be used to write the constraint in matrix form.The matrix   and vector   are used to represent the constraint as    ≤ 1   , where 1   is a column vector of   ones.Allow     to be the entry of   in the  ℎ row and  ℎ column.  can be expressed as The problem, including the group formulation can now be written as , (11)   Note that the group constraint only limits the number of selected incoming edges to be at most one.Due to the flow conservation constraint,  = , the number of chosen edges proceeding from each group will also be at most one.

Simple Example of the Group Constraint:
We continue with the previous graph example describing the group constraint depicted in Figure 3.There are two group constraints depicted in green and red in the figure.The first is around nodes  2 and  4 with the incoming edge set E 1 = { 1 ,  4 } and the index set I 1 = {1, 4}.The second is formed with nodes  3 and  4 with E 2 = { 2 ,  5 } and I 2 = {2, 5}.The group constraints can be written as

Data constraints
Nodes corresponding to an imaging or downlink task, will affect the data onboard the satellite.Let   represent the memory consumption at node  with   > 0 for imaging tasks,   < 0 for downlink tasks, and   = 0 for all non-task nodes.The data onboard spacecraft  at time  is denoted as    .The data vector over time for spacecraft  is denoted as   and all data vectors are combined in the vector , i.e., where   is the total number of time steps being considered,   is the number of spacecraft, and   is the total number of data variables.
The data at time  will depend upon the data at the previous time step and the edge taken to arrive at time .
The set    ⊂  represents the nodes associated with spacecraft  at time .As each node has an associated time value, two observations can be made about the structure of the graph: 1.The network flow constraint will restrict the number of edges entering    to be at most one.
2. The edge taken to arrive at time  will come from a node in    where  ≤ .
Denote the index set    as the set of edge indices corresponding to edges that terminate at a node in    .Recall that   denotes the index of the node where edge  terminates (i.e., if   = (  ,   ), then   = ).The updated memory consumed at time  can then be written as where  ,0 is the initial data storage.Note that multiplying   by    ensures that the contribution of node   to the total data is only received when an edge is chosen that enters node   .
The variables of optimization are now a combination of the edge variables and dynamic variables,  = [    ] ∈ R (  +  ) .Another index mapping,    , is used to map the data variables to their respective index in the optimization vector such that    =    .
The relation in (13) can then be written in matrix form for spacecraft  by defining the matrix where where The aggregate data update constraint can then be written as Note that each data element is assigned a zero utility (i.e.,   ,  = 0 ∀, ).The optimization problem with group and on-board memory constraints can then be written as

Example graph formulation
This section will now provide an example of the problem formulation for a single satellite performing imaging only (i.e., no downlink operations).A later example will address the coupled imaging and downlink scenario.The graph shown in Figure 4 will now be formulated using the network flow method outlined above.The edges are labeled with the letter 'x' and a numeric subscript to indicate the edge number.A numeric value is present with each edge to show the utility score associated with that edge.The grouping of nodes within each ellipse represent discrete time steps from  1 to  6 and show that only one node within each group can be visited within that time step since the satellites are modeled as satellites that can only perform one task at a given time.This application of grouping constraints is critical to the planning technique presented in this work and unique to the literature previously referenced.The resulting  matrix is shown in Figure 5.
Solving the problem, using the network flow implementation, results in the edge selection sequence of 2, 8, 15, 20, 23 as shown in Figure 4.The edges emphasized with green illustrate the optimal path through the graph when only considering imaging opportunities.
The network flow concept is now presented in the presence of imaging and downlink tasks within a memoryconstrained satellite environment.Consider the same graph from above but now with nodes 5, 8, and 9 be-  coming downlink opportunities instead of imaging operations.To properly create the   matrix, we must specify the amount of memory associated with each node, or edge, in the graph.
The memory associated with each operation is shown overlaid on the graph as purple numbers in Figure 7.
Note that memory associated with imaging operations are applied at the node but downlinking operations bookkeep memory on the edge between downlink nodes due to the time discretization applied in the problem.The actual amount of memory downlinked however, is the product of the length of the edge and the downlink rate.This ensures proper accounting for edge length since downlinking occurs over extended periods of time and imaging operations are modeled as instantaneous collections.Memory for downlinking is computed by multiplying the edge length (time) by the downlink rate of the spacecraft.Notice that the start and end nodes are assigned a memory value of 0 while the edges between nodes 5, 8, and 9 have all been allocated a memory amount of -11 to denote downlink to the ground terminal.The resulting  1  matrix (1 for a single satellite) is shown in Figure 6.
For this initial example, the memory limit is set to a value of 20.Note that this limit applies to every time step within the planning graph and ensures that the satellite does not collect more image data than it can save on-board.Solving this problem again with a memory limit set results in a new sequence through the graph nodes.This new sequence is shown in Figure 7. Note, that the new graph solution deviates from the original and passes through nodes 8 and 9 to reduce the amount of image data on-board the vehicle.The resulting data storage at each time step, from  1 to  6 , is given in Table 3.

GRAPH PRUNING
Recall that an edge is created between two nodes if the slew rate to achieve that transition is below a threshold.Even if two targets have very different slew angles, if they are located far enough apart in time then an edge will exist.The number of edges for the graph can quickly become very large with increased planning horizons or numbers of imaging targets due to the physics of satellite motion.In reality the only edges that are naturally disallowed, in Section 3.3, are those that transition between nodes that are located very close in time.While beneficial for expressing possible transitions, the complexity of the graph (and thus the number of decision variables and constraints considered) can quickly become unwieldy.
For example, consider the active target tasks depicted in Figure 8.Given 2,814 targets, 41 of which are visible to at least one satellite, the number of edges for three satellites over a planning horizon of 94.616 minutes is 7,440.Shown in the figure is a depiction of the resulting adjacency matrix with cyan depicting a possible edge.Since the goal is to realize a large number of target assignments, the density is not actually needed.Taking long edges results in skipping over a large number of possible imaging targets, making it an undesirable path.In fact, most optimal plans track closely to the diagonal, corresponding to satellites imaging as many targets as possible.An approach is now presented which can drastically reduce the graph complexity with little change in the resulting outcome.
A naive approach to reducing complexity is to eliminate all edges beyond a given maximum edge length.While it can significantly reduce the number of edges, it can also reduce the connectivity of the graph.A particular target may become disconnected from feasible paths from the start node to the end node.A group of nodes may be completely cutoff if the maximum edge length is too small.For example, lighting constraints and oceans can create large gaps in time between viable targets.Instead of solely using a maximum edge length, a solution using "backbone" nodes has been developed.Similar to the start and end nodes, the backbone node has no associated task and has a zero score.Such nodes are  placed between groups of nodes that would be disconnected due to the maximum edge length constraint.Nodes corresponding to each satellite are sorted in time and a linear search performed to find gaps larger than the maximum edge length.A backbone node is then placed in the temporal center of the gap and acts as a bridge between the separated groups.Edges are also added between the backbone nodes, when needed, to preserve connectivity.This method greatly reduces the number of edges, and thus the overall complexity of the graph.
Results for the backbone trimming approach using the example depicted in Figure 8 can be seen in Figures 9, 10, and 11.The horizontal axis in each figure shows the maximum edge length allowed.Each plot shows the results using an untrimmed graph, the naive trimming, and the backbone trimmed approach.Figure 9 shows that with a relatively small maximum edge length of 100, the same utility was achieved using the backbone node approach, but the naive trimming approach requires a length of roughly 3,000.The result is the ability to achieve nearly the same utility with a fraction of the computation time.Note that the number of edges, as depicted in Figure 10, has a strong correlation to the computation time.Also note that there is a small trade-off as the backbone approach does introduce additional nodes, becoming burdensome for very small maximum edge lengths.Figure 11 shows that the backbone approach achieves a similar reduction in graph complexity while maintaining a very similar plan to the untrimmed graph.

SIMULATION AND RESULTS
To fully examine the network flow constellation planning concept discussed throughout this paper, a realistic orbital flight scenario is analyzed in this section.The orbital scenario is first specified and then the results from planning are presented.A 500 km sun-synchronous orbit was chosen to provide realistic access periods to both targets and ground stations within the view of the Earth imaging satellites.The epoch for this scenario is August 1, 2021 at 18:00:00 UTC with accompanying orbital elements specified in Table 4.This orbit represents the seed satellite of a Walker delta constellation of 100 satellites and 25 orbital planes with no relative phasing between neighboring planes.A Walker delta pattern constellation provides a simple, parameterized method for conveying the orbital geometry for a constellation of many satellites.These parameters specify the inclination, total number of satellites, number of orbital planes, and the relative phasing between satellites in neighboring planes.This phasing angle is the angle difference, in the direction of orbital motion, from the ascending node to the closest satellite, when a satellite in the next westerly plane is at its ascending node.Again, in this example the relative phasing is 0°.A Walker delta pattern constellation starts with a seed satellite and generates the entire constellation from its orbital information.For more information about Walker constellations, see [38] .The ground terminals supporting downlinking operations are located in Alaska, Antarctica, Australia, Hawaii, and Norway.A total of 2,814 imaging targets were considered during planning, with these targets being approximately evenly spaced across the Earth' s land mass as depicted in Figure 12.

Results
The network flow constellation planning algorithm discussed throughout this paper was run on the mission scenario described previously to arrive at a constellation flight schedule.Several metrics were collected for each of the satellites to assess their performance and help evaluate the constellation planning routine.The metrics of particular interest in the DRM are satellite utility, number of targets collected, number of targets downlinked, and overall satellite memory state.These metrics are shown over the planning horizon timeline of approximately 3 hours (2 orbit revolutions) in Figure 13.The figure illustrates the performance that can be anticipated when using the network flow constellation planning routine presented.For context, Figure 14 also shows the satellite constellation with orbit trajectories, imaging events, and downlink operations illustrated.The results shown in Figure 13 provide insight into specific satellite assignments with red dashed lines representing individual satellites while the blue line in each plot highlights the collective performance across the constellation.The utility score, number of targets collected, and number of targets downlinked trend in an upward fashion since these metrics are cumulative and continue to increase as the satellites perform the imaging and downlinking operations over the simulation period.The satellite memory state plot is more variable as the satellites balance on-board memory capacity with new imaging collection opportunities and downlink operations.The memory plots show some satellites reaching a plateau in consumed on-board memory, indicating that imaging or downlink opportunities were either not selected or unavailable during that time, while other satellites show both imaging and downlinking operations occurring.The simulation was initiated with all satellites having all memory available on-board.Due to this, the global blue line trends upward for this simulation as memory across the constellation is consumed with imaging operations across the fleet.
Planning the imaging and downlink operations for 100 satellites over a multi-revolution period on a desktop PC is a challenging endeavor, yet the network flow technique proved to be capable.This implementation shows the ability to arrive at an optimal plan for a large constellation of 100 satellites while not exceeding the individual satellite' s resources.Furthermore, the plots show a relative global consistency for the constellation while also identifying unique trends and behaviors for each satellite.This consistency is due to the fairly homogeneous distribution of satellites relative to the imaging and downlink opportunities across the globe and lends confidence to the viability and consistency of the resulting solution.
The utility score and targets collected plots show a strong trend upward, meaning that imaging continues to progress throughout the planning period.Similarly, downlink is being optimized across the constellation, as illustrated by number of targets downlinked and onboard memory state plots, to ensure delivery of the data and prevent on-board resources from being exceeded.These results highlight the power behind the network flow method and support further research into this area.In particular, further research is warranted into additional methods for understanding the elapsed time from the collection of a specific image to the time that image information is delivered to the ground.The current implementation focuses on prioritized data collection and delivery but does not account for the value of quick delivery of image data.The latency metric is of particular interest to GEOINT missions due to the value in timely intelligence.Low-latency information allows planners to make time-critical decisions with greater confidence and is thus an important metric within the GEOINT community [39] .Global view of satellite planes and orbital configuration analyzed using 100 satellites in 25 planes with 5 ground terminals.In the left illustration, the red dots show where the satellite was located when an image was captured while the pink line illustrates the pointing at the time of collection.In the right illustration, the pink lines centered on the ground stations mark the downlink access periods.The targets colored in blue are those that were not considered during planning due to being in eclipse and not meeting the specified lighting constraints.

CONCLUSION
This paper has provided a detailed explanation of the constellation scheduling problem, within the context of the geospatial intelligence design reference mission, as well as a comprehensive overview of the graph method and network flow approach to arrive at a globally-optimal solution within the discretized problem space.A realistic constellation planning scenario was presented for 100 satellites, prosecuting over 2,800 targets, while accessing up to 5 ground station terminals.The results show tractability and trending behaviors of increasing utility by collecting image data and downlinking it to the ground, while keeping within the onboard resources of each satellite within the constellation.This paper has also shown the dramatic improvement in the time required to solve the geospatial intelligence constellation planning problem once graph pruning has been applied.This improvement is the first step in making the algorithm run fast enough to support feasible solution generation in an operational environment.Furthermore, the presented methodology illustrates how a variety of specific spacecraft operational constraints can be created using network flow theory.The problem formulation and solving methods described in this paper have far-reaching implications and provide a robust strategy for planning coordinated operations across a fleet of cooperative systems in a time-extended mission context.Future research will leverage the network flow capabilities and extend the problem space by introducing cross-schedule dependencies for direct interaction between satellites and support crosslinking of data between these space systems.This capability will allow satellites to share data between them, thus freeing up onboard resources, while also potentially improving the global utility performance of the constellation, as a whole.

Figure 1 .
Figure 1.The constellation planning process from mission configuration through command generation.

Figure 2 .
Figure 2. Depiction of a valid access window using elevation and azimuth constraints on the left and a discretized satellite trajectory through the window on the right.

Figure 3 .
Figure 3.A simple graph is presented with numbered circles as nodes and labeled edges.The right image illustrates two group constraints, one around nodes 2 and 4 and one around nodes 3 and 4.

Figure 4 .
Figure 4.The optimal path through the graph when solving this problem using the network flow approach outlined.The path through the nodes is given as edges 2, 8, 15, 20, and 22.

Figure 5 .
Figure 5.The resulting  matrix from the example graph.

Figure 6 .
Figure 6.The  1   matrix representing the flow of memory through the graph.

Figure 7 .
Figure 7.The path solution resulting when constraining the on-board memory to a value of 20.The path of edges is given as 2, 8, 14, 16, and 21.

Figure 8 .
Figure 8.On the left is shown the imaging target deck, the middle shows the adjacency matrix assuming a single satellite, and the right image shows the adjacency graph for three satellites.The adjacency matrix is visualized using a cyan pixel for a 1 and white pixel for a 0. The nodes are sorted first by time and then by satellite.The red lines show the optimal path through the graph.

Figure 9 .
Figure 9.The resulting utility score and computation cost from the scenario in Figure 8 using three satellites.The left image shows the edge length vs utility for the untrimmed, naively trimmed, and backbone trimmed graphs.The middle image is a zoomed in view of the left image.The right image shows the resulting computation time for the trimming.

Figure 10 .
Figure 10.The number of edges and nodes vs the trim length.

Figure 11 .
Figure 11.The result of trimming on the adjacency matrix of the graph.The left matrix shows the untrimmed graph.The middle matrix shows the trimmed graph without backbone nodes.The right matrix shows the trimmed graph with backbone nodes.The red lines on the plot show the planned path for each of the three satellites in this example.The max edge length for these graphs is set to 300.

Figure 12 .
Figure 12.Imaging target distribution across the Earth's land mass area with targets accessible during the planning period marked in green.Targets marked in red do not meet the lighting constraints due to being in eclipse during the planning period.

Figure 13 .Figure 14 .
Figure 13.These figures show the results from a simulation with 100 satellites in 25 planes with 2814 ground targets spread over the surface of the earth and 5 ground stations in Alaska, Antarctica, Australia, Hawaii, and Norway.The left axis (in blue) shows the collective satellite values and the right axis (in red) shows the values for individual satellites

Table 1 . Notation used for Planning Name Description Name Description
Incidence matrixIncidence matrix for satellite

𝑗
The   ℎ row and   ℎ column of  The   ℎ row and   ℎ column of   The   ℎ row and   ℎ column of The ℎ row of     ℎ element of RHS of flow constraint  RHS of flow constraint with one source, one sink, and unit flow Group constraint parameters V Nodes in group ,  ⊂  E Set of edges terminating at node in V I Set of edge indices in E  Group constraint matrix  Marginal memory requirement for node   Nodes associated with time   Edge indices that terminate at a node in   Aggregate RHS vector for the data constraint Indices and index mappings  Index of the node where edge  terminates  Mapping of  to index within  , ,

Table 2 . Definition of task constraints Value Description
lat/lon/alt Latitude, longitude, and altitude of the target min/max azimuth