Nomenclature
=  onorbit depot spare capacity, in units of modules 
=  spare demand, in units of modules per hour 
=  lead time, in units of hours 
=  number of customer satellite modules, in units of modules 
=  service time, in units of hours 
=  inbound travel time, in units of hours 
=  outbound travel time, in units of hours 
=  repair time, in units of hours 
=  extra service time due to stockout, in units of hours 
=  launch interval, in units of hours 
=  time until stockout, in units of hours 
=  waiting time until service completion, in units of hours 
=  waiting time in queue, in units of hours 
=  mean individual module failure rate, in units of failures per hour 
=  mean launch rate, in units of launches per hour 
=  mean systemlevel spare demand rate, in units of modules per hour 
=  fill rate 
=  fill rate requirement 
=  probability density function 
=  expected value 
=  LaplaceStieltjes transform of a function 
1 Introduction
Nowadays, there has been increasing interest in developing onorbit infrastructure systems that enable sustainable space exploration. Over the last several decades, research and development in autonomous and robotics systems have significantly raised the technology readiness level of robotic onorbit servicing (OOS) [NASA, JAXA, DARPA]. Engineers envision spacebased servicing infrastructures to provide refueling and repair services or to manufacture large structures directly in orbit. The recent trend of satellite modularization is also enabling “servicingfriendly” satellites that are composed of multiple small standardized structural modules [Sullivan2001, Hill2013, Kerzhner2013, Jaeger2013]. These OOS infrastructure systems and serviceable satellite designs are expected to be a gamechanging alternative to the current satellite industry [Saleh2002, Lamassoure2002, Saleh2003, Long2007, Saleh2001].
Traditionally, most OOS concepts studied in the literature have assumed a dedicated robotic spacecraft to repair or refuel the customer satellites [Yao2013, Zhao2016, Verstraete2018]. The servicer would visit a predefined set of satellites and be discarded in a graveyard orbit once the mission is over. The advantage of this concept is that the flight path of the servicer can be optimized before launch to provide the best value to the servicing operation. However, although such a “disposable" servicer is a realistic concept in the near term, it is not a sustainable solution for OOS in the longer term.
Alternatively, a more sustainable concept is to use a permanent reusable servicing infrastructure that responds to random failures [Long2007, SartonduJonchay2017]. One example of such a concept would be as follows: when a module fails, the OOS system dispatches a servicer to the satellite to replace the failed module with a new spare. In this case, a permanent OOS infrastructure would include a servicer replacing defective modules with new module spares, an onorbit depot storing the spare modules, and a series of launch vehicles supplying new modules from Earth to the depot on a regular basis. This concept can provide timely and responsive services to the random failures spontaneously and thus enable sustainable space exploration.
However, despite the potential advantage of a permanent responsive OOS infrastructure system, its design and operations planning becomes increasingly complex and challenging compared with the traditional dedicated servicer concepts. This is because the analysis of a responsive OOS system would need to consider the interactions between the infrastructure elements as well as the full supply chain of the spare modules, including the queue of the services and the inventory of the spare depot. We need an efficient model to analyze the design of the OOS system.
In response to this background, we develop a novel semianalytical model to evaluate the OOS system design and analysis. Given a set of distributed modular satellites (i.e., customer satellites) with random failure rates, the proposed model analyzes the performance of the servicing system considering the spares supply chain. Specifically, our goal is to analyze the mean waiting time before service completion for a given customer satellite failure and its tradeoff relationship with the depot capacity.
The developed semianalytical model leverages the queueing theory and inventory management methods. Queueing theory is the mathematical theory of waiting in lines; it models realworld queueing systems using distributions of customer arrival, service time, and queue discipline [LarsonOdoni]. Inventory management refers to the process of ordering, storing, and using inventory, such as raw materials, components, and end products; in this paper, inventory control methods are applied to find an optimal order policy that aims to balance inventory cost and demand shortage [SimchiLevi]. Thus, the developed model contains a set of coupled submodels: (1) a queueing submodel that models the mean waiting time for the module failures; and (2) an inventory control submodel that models the replenishment of the onorbit depot. The results from the semianalytical model are compared with simulation results for different realworld cases, and the accuracy of the proposed model is demonstrated.
A few remarks need to be made about the value and contribution of the semianalytical model developed in this paper. Conventionally, the analysis of space systems with random failures and repairs (e.g., OOS systems) has been performed using computationally costly simulations [SartonduJonchay2017, SartonduJonchay2017a, Sears2018]. However, as space systems become complex, designing and evaluating their performance using only simulations becomes computationally challenging. While simulations are effective for detailed design, there is a growing need for the development of more efficient, yet rigorous, analysis methods to enable quick performance evaluation for systems design and trade space exploration. Although analytical or semianalytical models have been recently introduced gradually in space systems design [Jakob, Zchen], this research direction remains largely unexplored. This paper introduces the first semianalytical approach with an integrated queueing and inventory model for complex space systems analysis. The model developed in this paper is expected to be a critical step in pushing forward this research frontier of analytical/semianalytical model development for space systems design.
The rest of the paper is organized as follows. Section 2 discusses the overview of the considered OOS architecture. Section 3 explains the main contribution of the work: a semianalytical model for OOS architecture. Section 4 assesses the results from the proposed model using simulations with a realistic application example, and Section 5 concludes the paper.
2 Overview of OnOrbit Servicing Infrastructure
This paper considers an OOS architecture to provide modulereplacement services to a set of customer satellites distributed in space (e.g., at the geosynchronous equatorial orbit). The customer satellites are modular, where each module can fail and be replaced independently. In this study, we consider all modules to be identical for simplicity, but the developed model can be extended to multiple types of modules. The servicing architecture is comprised of three main components: the servicer, the onorbit depot, and the launch vehicles. The onorbit depot stores modules which are brought to customers by the servicer when failures (i.e., demands) happen. The depot is kept at the same orbit as the customers and is refilled with modules by a series of launch vehicles. The servicing architecture can be seen in Fig. 1.
The overview of the concept of operations is as follows. The servicer remains docked to the depot while awaiting a customer satellite module failure. Once a customer satellite module fails, the servicer performs a phasing maneuver to the failed module. The servicer performs the repair operation, defined as the replacement of a failed module with a spare module, and then returns to the depot. For multiple customer failures, the services are performed on a firstcomefirstserved (FCFS) basis; namely, customers receive service in the order in which they failed, regardless of their relative orbital positions to the servicer. The servicer can only carry one spare module, and thus it must return to the depot to resupply and refuel itself (these operations are considered to be instantaneous) prior to the next repair trip. The spare inventory in the depot is monitored regularly, and the launch vehicle visits the depot for its resupply following a stochastic launch schedule.
3 SemiAnalytical Model
This section provides an overview of the problem statement, the details of the developed model, and the proposed formulations based on the model.
3.1 Problem Statement
We develop a novel semianalytical model to evaluate the OOS system design and analysis without relying on computationally costly simulations. As the first step, this subsection converts the concept of operations of the considered OOS system into a mathematical problem.
The overall concept of operations can be expressed in a block diagram in Fig. 2. Here, the queue for the services needs to be analyzed rigorously to capture both the service time and the waiting time in the queue before the service operation starts. We define the service time and the waiting time (i.e., the time between a module failure and its repair completion) as follows:
(1) 
(2) 
where is the delay of the service when the depot is found to be out of stock; and are the outbound and inbound (i.e., return) transportation times, which depend on the relative positions between the customer satellites and the depot; is the repair time, which is assumed to be a fixed value in this paper but can be varied if needed; and
denotes the waiting time in the queue. We can find the probability distribution of
and from the (known) distribution of the positions of the customer satellites. is derived from the inventory control method, and is derived from the queueing method. Note that, technically, is not a service time, but we consider it as part of the service time so that the queueing method can naturally consider inventory management with little modification.The considered inventory control strategy of the depot can be modeled as a modified orderupto policy. Namely, over every exponentially distributed time interval
, we order units from the ground, where is the capacity of the depot, is the current inventory level, is the replenishment on the way (i.e., orders being processed). This exponentially distributed time interval corresponds to the launch frequency of the rocket with a mean launch rate , which is known to be a good approximation [SartonduJonchay2017]. A constant lead time will be added to account for the processing of the order, manufacturing of the units, and loading of the units to the rocket. Thus, time steps before each launch opportunity, the inventory level is checked and the orders start to be processed. Unmet demand due to stockout is backordered and delivered in the following interval. Fig. 3 summarizes the considered inventory control policy.A key consideration for the inventory performance is the fill rate , representing the percentage of the spare demand that is met without stockout. We assume that we are given the fill rate requirement (typically close to one), and our goal is to analyze the mean waiting time and the depot capacity . The tradeoff between and corresponds to a costperformance tradeoff. We prefer a lowcost and highperformance OOS system, which can be interpreted as a system with a small and a small ; however, a smaller indicates a larger , and thus the tradeoff between these two needs to be considered. Particularly, we are interested in the behavior of vs. around the point of diminishing returns (i.e., knee point), which represents the depot capacity above which little additional saving is expected in the mean waiting time . This region corresponds to where is close to one, which is true for a realistic reliable OOS system. This concept is shown in Fig. 4.
3.2 Proposed Model
This subsection introduces the proposed model to analyze the performance of the OOS system. The developed model leverages the queueing theory and inventory control methods; the integrated queueing and inventory model is summarized Fig. 5. The inputs to the model include the probability distributions of the transportation and repair times (), the individual module failure rate , the number of modules , the length of lead time , the mean launch rate , and the required fill rate . The integrated queueing and inventory control model takes those inputs and derives the mean waiting time and the depot capacity . The queueing and the inventory control submodels interact as follows: the queueing submodel generates the spare demand (i.e., module failure) distribution , which is fed into the inventory control submodel to calculate the probability distribution of the extra servicing time due to stockout ; this distribution is fed back into the queueing submodel. Since it is difficult to analytically handle the spare demand distribution due to its statedependent nature, we approximate this by a Poisson process with the corresponding mean demand rate . In this way, can be characterized using only , which is part of the output from the queueing submodel. (Note that this value is not just the consolidation of the individual module failure rate because the system failure rate is statedependent; see Section 3.2.1.) This approximation is demonstrated to perform well in later numerical simulations.
The following subsubsections introduce the details of the queueing submodel and inventory control submodel in the integrated model.
3.2.1 Queueing Submodel
The queueing part of the OOS problem can be modeled using a finitesource queue. A finitesource queue represents the case where the number of customer modules is finite. Consequently, as the number of failed modules increases, the number of active modules decreases, thus decreasing the systemlevel failure rate. This fact makes the failure rate statedependent (i.e., the failure rate depends on the number of active modules), and thus makes the problem challenging.
Using Kendall’s notation, the considered queue is written as . The meaning of each letter is as follows:

: The arrival process is Markovian.

: The service time distribution is general.

: The number of servicers is one.

: The number of failures allowed to in the queue is .

: The size of the source population is .
The general solution for the queue can be found using the LaplaceStieltjes transform of the service time distribution. Denoting the LaplaceStieltjes transform of a function as , the mean spare demand rate becomes:
(3) 
and the mean waiting time is:
(4) 
where
(5) 
(6) 
Derivation of the above equations and its computationally efficient algorithms can be found in Refs. [Takacs1962] and [Gupta1994]. Here, the mean and the probability distribution of the service time, and , can be found using Eq. (1), where the distributions of , , and are known and the distribution of can be found using the inventory analysis (see Eq. (12)). Also, in Eq. (3) is used to generate the demand distribution with a Poisson assumption as discussed above:
(7) 
This is then fed into the inventory control analysis.
Note that an implicit assumption for this queueing analysis is that, since the servicer’s return transportation time is included as part of the service time, the replaced module does not resume normal operation (i.e., module failure process does not restart) until the servicer returns to the depot after its service. This approximation is reasonable in practice because the modules’ mean time between failures (MTBF) is much longer than the transportation time.
3.2.2 Inventory Control Submodel
The inventory control analysis for the OOS problem contains two parts: finding given and finding the . In addition to the Poisson assumption for , we make two assumptions for the inventory model:

When the depot becomes stockout, the backordered units are delivered together with the capacity of units. This may not apply in reality due to the rocket capacity constraint; however, it is a reasonable approximation when the solutions are near the knee point, where a stockout happens rarely. ^{6}^{6}6In later simulations, we consider realistic cases where we only order up to units even when the inventory is zero; nevertheless, the proposed model is demonstrated to approximate the simulation results well. This assumption prevents the queue from blowing up in the developed model.

The expected backorder and its impact on the servicing time are computed assuming no resupply during the interval and the lead time. This is a conservative assumption.
For the orderupto policy, following the inventory control literature [Cachon], the fill rate is defined as follows:
(8) 
where the exponential launch interval distribution is as follows:
(9) 
and the demand distribution can be found in Eq. (7). The numerator of the second term in Eq. (8) is the expected backorder over a launch interval and a lead time, whereas the denominator is the expected demand over a launch interval. The reason why the numerator considers both the launch interval and the lead time (instead of only the launch interval) can be intuitively explained with the following example. Consider a case where an order is made at a point where the inventory level is units. In this case, the order would be delivered after time steps where the inventory level is units, where is the additional units of demand during the lead time . Because the order only delivers units using the information when the order was made, the actual inventory right after the replenishment delivery is instead of the full capacity . Therefore, a stockout (and thus a backorder) happens when the demand over the next launch interval, denoted as , exceeds ; this corresponds to when the summation of the demand over an interval () and the demand over the lead time () exceeds . This explains why the fill rate computation needs to take into consideration the backorder over both the launch interval and the lead time. For a more rigorous derivation of the fill rate expression in Eq. (8), see Ref. [Cachon].
A closedform expression for (Eq. (8)) using a given demand rate can be derived and is summarized in Appendix A. This expression allows us to evaluate the fill rate given a depot capacity . With this expression, we can find the minimum that satisfies a given fill rate requirement :
(10)  
s.t.  (11) 
where is given by Eq. (8). Since is a monotonically increasing function of , the solution to this problem can be found simply by iteratively incrementing until is satisfied.
Using the depot capacity , we can derive the expression for the additional service time due to stockout , which is then fed back to the queueing analysis. corresponds to the extra "service" time that the first repair that finds the depot to be out of stock needs to wait additionally before starting the operation. As long as it takes the same or a longer time to experience demands () than the sum of launch interval and the lead time (), there is no stockout. However, when that is not the case (), then the first demand that finds the depot to be out of stock experiences extra service time. (Note that this is a conservative approximation; in reality, the extra service time would only be added after all failures in the queue are serviced.) Only the first one out of demands on average per launch interval will be affected; for all the remaining demands, no extra service time will be considered. Thus, can be written as follows:
(12) 
where and
themselves are also random variables following the probability density functions
and , respectively. Note that the Poisson demand assumption makes the time until failures (i.e., ) represented with an Erlang distribution:(13) 
A closedform expression for (Eqs. (3)(6)) using this (Eq. (12)) and a given demand rate can be derived and are summarized in Appendix B. Note that this expression for can lead to an expression for , which is used in this very expression for ; thus, this system of equations needs to be solved concurrently.
3.3 Proposed Formulations
With the proposed queueing and inventory control submodels, we develop two integrated formulations for the OOS system. The first one focuses on approximating the knee point computationally efficiently given a high fill rate; the second one aims to approximate a series of solutions around the knee point more precisely.
3.3.1 Formulation 1: Finding the Knee Point
The goal of Formulation 1 is to computationally efficiently approximate the knee point given a high fill rate. To this end, we decouple the loop between queueing analysis and the inventory control analysis by applying the following approximation:
(14) 
Thus, is not needed for the queueing analysis. With this approximation in Eq. (14), we can separately use the queueing analysis to find the mean waiting time and mean demand rate via Eqs. (3)(4), and then use the inventory control analysis to find the depot capacity to satisfy the fill rate requirement via Eqs. (10)(11); no feedback loop is needed. The computed mean waiting time will be the case with no stockout, and thus an optimistic bound.
3.3.2 Formulation 2: Finding the Mean Waiting Time vs. Depot Capacity Relationship Around the Knee Point
The goal of Formulation 2 is to find a series of solutions around the knee points over different . In this case, we consider the integrated model and solve the coupled queueing equations (Eqs. (3)(4)) and inventory control equations (Eqs. (10)(11)) concurrently; the former equations need from the latter, whereas the latter equations need (i.e., in the form of ) from the former (see Fig. 5). These coupled equations generally cannot be solved fully analytically; a standard numerical solver such as the fsolve
function in MATLAB can be leveraged.
The solution from Formulation 1 (which can be efficiently obtained) can be used as an initial value for the fsolve
function. This formulation enables us to find both and over different around the knee point.
4 Application Example
This section applies the proposed method to an example case with realistic parameters and assesses the accuracy of the proposed model with simulations. The considered example case contains 10 satellites with 5 modules each that are evenly distributed over the geosynchronous orbit. The parameters are listed in Table 3. Three cases for the module MTBF are considered.
Parameter  Value  

Repair time  4 hours  
Launch lead time  2160 hours  
Mean launch interval  1213.4 hours  
MBTF per module  20000, 10000, 2000 hours 
For the considered application case, the developed semianalytical model is used to derive: (1) the estimated knee point (i.e., Formulation 1); (2) the mean waiting time vs. depot capacity curve (i.e., Formulation 2).
^{7}^{7}7For this illustrative example, fsolve function in MATLAB is used with the function tolerance of , the optimality tolerance of , and the step tolerance’ of for Formulation 2. To assess the accuracy of the proposed semianalytical model, we use the depot capacity derived from the formulations and evaluate the corresponding waiting time of the system via simulations; the mean waiting time results from the simulations are compared with those from the semianalytical model. For each given set of the MTBF and the depot capacity, 50 simulation runs are performed over a time horizon of 200000 hours (i.e., 22 years). Note that both formulations can provide the same value of the depot capacity for a given set of the MTBF and the fill rate requirement, in which case the same set of simulation results is used for the evaluation of both formulations. The simulation methods used in this paper are based on Ref. [SartonduJonchay2017]. We evaluate the cases with fill rate requirements , expecting a good performance of the proposed formulations when is close to one.The results of the semianalytical model and the simulations are shown in Figs. 6 and 7. For the simulation results, the averaged waiting times over the time horizon for each 200000hour run are evaluated for the corresponding depot capacities and are shown as individual black dots in Figs. 6 and 7; additionally, their means over the 50 simulations for each depot capacity are connected by a black solid line. For the semianalytical results from Formulation 1, we derive the mean waiting time with no stockout (i.e., infinite capacity) and the depot capacity given the fill rate requirements. These results are shown as the crosses for each given fill rate as well as the dashed line connecting them in Fig. 6. For the semianalytical results from Formulation 2, we derive the curve of mean waiting time vs. depot capacity near the knee points; these results are shown as the crosses for each fill rate as well as the dashed line connecting them in Fig. 7.
The quantitative comparison of the results from the model and simulations is also shown in Table 4 for Formulation 1 and Table 5 for Formulation 2. Note that the error between the model and the simulations is caused by both the approximation made in the model and the randomness in the simulations due to the finite number of runs.
MTBF, hours  , hours  , hours  Error in  

20000  0.8  12  30.5  387.8  92.1% 
0.85  13  30.5  288.5  89.4%  
0.9  15  30.5  145.7  79.0%  
0.95  17  30.5  90.2  66.1%  
0.99  23  30.5  39.9  24.5%  
0.995  25  30.5  37.3  18.2%  
0.999  31  30.5  31.0  1.4%  
0.9995  33  30.5  31.7  3.6%  
10000  0.8  23  35.5  367.0  90.3% 
0.85  24  35.5  293.4  87.9%  
0.9  27  35.5  176.8  79.9%  
0.95  32  35.5  92.6  61.7%  
0.99  42  35.5  48.4  26.8%  
0.995  47  35.5  40.2  11.8%  
0.999  57  35.5  36.9  3.9%  
0.9995  62  35.5  36.7  3.3%  
4000  0.8  52  65.8  464.7  85.8% 
0.85  57  65.8  334.4  80.3%  
0.9  63  65.8  216.9  69.7%  
0.95  74  65.8  162.1  59.4%  
0.99  98  65.8  85.0  22.5%  
0.995  109  65.8  70.3  6.4%  
0.999  134  65.8  67.8  2.9%  
0.9995  144  65.8  67.5  2.5% 
MTBF, hours  , hours  , hours  Error in  

20000  0.8  12  306.1  387.8  21.1% 
0.85  13  232.1  288.5  19.5%  
0.9  15  140.5  145.7  3.6%  
0.95  17  91.5  90.2  1.5%  
0.99  23  41.3  39.9  3.4%  
0.995  25  36.6  37.3  1.9%  
0.999  31  31.6  31.0  2.1%  
0.9995  33  31.2  31.7  1.7%  
10000  0.8  22  330.1  437.8  24.6% 
0.85  24  245.7  293.4  16.3%  
0.9  27  171.3  176.8  3.1%  
0.95  32  96.6  92.6  4.3%  
0.99  42  48.4  48.4  0.1%  
0.995  47  41.4  40.2  3.1%  
0.999  57  36.8  36.9  0.4%  
0.9995  62  36.1  36.7  1.7%  
4000  0.8  48  449.4  588.6  23.6% 
0.85  54  323.5  395.4  18.2%  
0.9  61  233.6  251.4  7.1%  
0.95  73  143.5  141.6  1.3%  
0.99  98  81.6  85.0  3.9%  
0.995  109  73.6  70.3  4.7%  
0.999  134  65.6  67.8  3.3%  
0.9995  144  66.7  67.5  1.2% 
From the results, we have the following findings for the developed formulations:

Formulation 1 evaluates the bestcase scenario when there is no stockout (i.e., infinite depot capacity). For the tested cases, its mean waiting time evaluation achieves the accuracy of when . This formulation is beneficial when we design a highly reliable OOS system or are interested in evaluating the ideal scenario.

Formulation 2 explores the relationship between the mean waiting time and the depot capacity near the knee point. It achieves the accuracy of when . This formulation is beneficial when we explore a tradeoff between performance and cost.
Note that, in practice, the OOS system is designed to have a high order fill rate so that the services can be completed within reasonable waiting times; therefore, the approximation in the developed formulations serves the purpose. The above results demonstrate the accuracy and utility of the proposed semianalytical model for the considered OOS application.
One benefit of the proposed model is that it does not require computationally costly simulations for evaluation. For the tested cases, running a set of 50 simulation runs over a 20000hour time horizon for each case takes more than 30 hours with Python 3.6, whereas the semianalytical model completes its evaluation within 0.1 seconds for Formulation 1 and within 1.5 minutes on average for the cases tested for Formulation 2 (7.5 minutes for the worst case tested) with MATLAB R2019a, all on an Intel Core i78650U CPU @ 1.90GHz platform with 16GB RAM. Although the exact computational time may vary depending on the implementation, the computational cost saving by the developed model is evident. The developed model provides an efficient highlevel design analysis and optimization method with sufficient accuracy for earlystage systems design, complementing the existing costly simulation techniques for detailed design.
5 Conclusion
This paper develops a novel semianalytical model for OOS system analysis based on queueing theory and inventory management methods. We consider an OOS system that provides responsive services to the random failures of a number of modular customer satellites in orbit. The considered OOS architecture is comprised of a servicer that provides modulereplacement services to the customer satellites, an onorbit depot that stores the spare modules, and a series of launch vehicles to refill the depot. The developed model is capable of analyzing the queue of the service operations as well as the logistics of the spares, evaluating the mean waiting time before service completion for a given failure and its relationship with the depot capacity. The case study shows that the results from the semianalytical model approximate the solution well without computationally costly simulations.
Although this paper uses a simple OOS application case to demonstrate the value of the developed model, the developed general model can be applied to more complex OOS problems with reasonable modifications. Possible extensions include: (1) employment of a servicer that can carry multiple spares (i.e., servicing multiple customers in one trip); (2) consideration of multiple types of spares, repair operations, or even servicer spacecraft; and (3) consideration of safety stock at the depot to minimize the stockout risk. We expect that this paper opens up a broad range of applications of semianalytical models to OOS system design.
Acknowledgment
The authors thank Brian Hardy for running the simulations needed for this paper and Hang Woon Lee for reviewing this paper.
Funding Sources
This work is partially supported by the Defense Advanced Research Projects Agency (DARPA) Young Faculty Award D19AP00127. The content of this paper does not necessarily reflect the position or the policy of the U.S. Government, and no official endorsement should be inferred. Approved for public release; distribution is unlimited.
Appendix A: ClosedForm Expression for
A closedform expression for the fill rate in Eq. (8) can be derived as follows:
Appendix B: ClosedForm Expression for
This appendix summarizes a closedform solution for the mean waiting time in Eqs. (3)(6). The only complication is obtaining the analytical expression for and . can be expressed as:
(15) 
Similarly, leveraging the convolution theorem for the LaplaceStieltjes transform, can be expressed as:
(16) 
The terms for , , and in Eqs. (15)(16) can be derived using orbital mechanics; deriving this analytical expression is trivial and thus omitted because we know the exact locations of the customer satellites. Note that this expression is already sufficient for Formulation 1.
The term for in Eqs. (15)(16), needed for Formulation 2, is more challenging. We first derive and use that to find .
(17) 
Focusing on ,
Thus, we have obtained a closedform expression for .
Using , can be found as follows:
Thus, we have obtained all terms needed for a closedform expression for .