Orbital Facility Location Problem for Satellite Constellation Servicing Depots

TL;DR

This paper adapts the Facility Location Problem to optimize the placement of on-orbit servicing depots for satellite constellations in high-altitude orbits, considering transfer costs, launch costs, and orbital dynamics.

Contribution

It introduces an orbital facility location model using binary linear programming and Lyapunov feedback control for transfer calculations, tailored for high-altitude satellite servicing.

Findings

Optimized depot placement for Galileo and GPS constellations.

Demonstrated cost-effective on-orbit servicing architectures.

Validated approach with real satellite constellation data.

Abstract

This work proposes an adaptation of the Facility Location Problem for the optimal placement of on-orbit servicing depots for satellite constellations in high-altitude orbit. The high-altitude regime, such as Medium Earth Orbit (MEO), is a unique dynamical environment where existing low-thrust propulsion systems can provide the necessary thrust to conduct plane-change maneuvers between the various orbital planes of the constellation. As such, on-orbit servicing architectures involving servicer spacecraft that conduct round-trips between servicing depots and the client satellites of the constellation may be conceived. To this end, orbital facility location problem is a binary linear program, where the costs of operating and allocating the facility(ies) to satellites are considered in terms of the sum of Equivalent Mass to Low Earth Orbit (EMLEO), is proposed. The low-thrust transfers between the facilities and the clients are computed using a parallel implementation of a Lyapunov feedback controller. The total launch cost of the depot along with its servicers, propellant, and payload are taken into account as the cost to establish a given depot. The proposed approach is applied to designing on-orbit servicing depots for the Galileo and the GPS constellations.

Tables

$a$  =	semimajor axis
${𝒄}_{ij}$  =	facility allocation cost
$D$  =	Demand of servicing trips
$e$  =	eccentricity
${𝒇}_j$  =	facility usage cost
$g_0$  =	gravitational acceleration
$i$  =	inclination
$i_{\text{sp}}$  =	specific impulse
$m_{d,\text{dry}}$  =	Depot dry mass
$m_{l,\max }$  =	Maximum launch mass of launch vehicle
$m_{s,\text{dry}}$  =	Servicer dry mass
$m_{s,L}$  =	Servicer payload mass
$r_a$  =	apogee radius
$r_p$  =	perigee radius
$𝒙$  =	Spacecraft orbital elements
$X_{ij}$  =	binary allocation variable
$Y_j$  =	binary usage variable
$Z$  =	initial-to-final mass-ratio
$\mu$  =	gravitational parameter
$\Omega$  =	right-ascension of ascending node
$\omega$  =	argument of periapsis

Table 2. Table 2 : OFLP and Q-Law parameters for Numerical Experiments

Parameter	Value
Maximum launch mass $m_{l,\max }$, $\mathrm{kg}$	12,950
Launch vehicle parking orbit radius $r_0$, $\mathrm{km}$	6,578
Launch vehicle $I_{\mathrm{sp},l}$, $\sec$	457
Payload mass $m_{s,L}$, $\mathrm{kg}$	100
Number of trips $D$	1 / 2
Servicer dry mass $m_{s,\text{dry}}$, $\mathrm{kg}$	500 / 1,000
Servicer $I_{\mathrm{sp},s}$, $\sec$	1,790
Servicer thrust, $N$	1.74
Depot dry mass $m_{d,\text{dry}}$, $\mathrm{kg}$	1,500 / 2,000 / 2,500
Depot $I_{\mathrm{sp},d}$, $\sec$	320
Maximum transfer time, $\mathrm{days}$	300
Minimum transfer radius $r_{p\min }$, $\mathrm{km}$	6,878
$W_p$	1
$W_{\text{œ}}$	[1,1,1,1,1]
$\sigma$	3
$\nu$	4
$\zeta$	2
$k_{r_p}$	1

Table 3. Table 3 : Candidate depot slots for GPS and Galileo Constellations

Orbital element	Value (min:increment:max)	Number of slots
Semimajor axis $a$, $\mathrm{DU}$ ($1\mathrm{DU}=26,560\mathrm{km}$)	0.3:0.05:1.1	17
Eccentricity $e$	0:0.05:0.6	13
Inclination $i$, deg	50:1:58	9
RAAN $\Omega$, deg	0:30:330	12
Argument of perigee $\omega$, deg	0	1
Total	-	23868

Table 4. Table 4 : CPLEX performance solving OFLP for GPS & Galileo servicing facilities

		$D=1$				$D=2$
$m_{s,\text{dry}}$, $\mathrm{kg}$	$m_{d,\text{dry}}$, $\mathrm{kg}$	Status	Gap	Solve time, $\sec$		Status	Gap	Solve time, $\sec$
500	1,500	OPTIMAL	0.0	655.929		OPTIMAL	0.0	174.511
	2,000	OPTIMAL	0.0	155.637		OPTIMAL	0.0	119.195
	2,500	OPTIMAL	0.0	156.419		OPTIMAL	0.0	133.745
1,000	1,500	OPTIMAL	0.0	295.819		OPTIMAL	0.0	257.044
	2,000	OPTIMAL	0.0	461.123		OPTIMAL	0.0	281.117
	2,500	OPTIMAL	0.0	787.712		OPTIMAL	0.0	2,709.114

Table 5. Table 5 : Summary of optimal GPS & Galileo servicing depots from OFLP instances with m s , dry = 500 subscript 𝑚 𝑠 dry 500 m_{s,\text{dry}}=500   kg times absent kg \text{\,}\mathrm{k}\mathrm{g}

		$D=1$								$D=2$
$m_{d,\text{dry}}$, $\mathrm{kg}$	No.	$a$, $\mathrm{DU}$	$e$	$i$, $\mathrm{deg}$	$\Omega$, $\mathrm{deg}$	${𝒦}_a$	$m_{d,\text{wet}}$, $\mathrm{kg}$	Depot EMLEO, $\mathrm{kg}$	Total EMLEO, $\mathrm{kg}$		$a$ $\mathrm{DU}$	$e$	$i$ $\mathrm{deg}$	$\Omega$ $\mathrm{deg}$	${𝒦}_a$	$m_{d,\text{wet}}$ $\mathrm{kg}$	Depot EMLEO, $\mathrm{kg}$	Total EMLEO, $\mathrm{kg}$
1,500	1	0.80	0.20	56	30	16	6,015	9,470	38,079		0.95	0.10	55	30	16	9,612	15,416	59,833
	2	0.60	0.55	57	90	5	2,758	4,312			0.60	0.55	57	90	5	3,973	6,211
	3	0.60	0.55	55	150	14	4,845	7,575			0.60	0.55	55	150	14	8,147	12,737
	4	0.55	0.50	54	210	4	2,501	3,799			0.60	0.55	53	210	4	3,343	5,227
	5	0.60	0.55	57	270	13	4,682	7,320			0.60	0.55	57	270	13	7,821	12,226
	6	0.55	0.50	54	330	7	3,688	5,603			0.95	0.05	55	330	7	5,048	8,016
2,000	1	0.75	0.25	55	30	16	6,734	10,546	42,294		0.90	0.10	55	30	16	10,475	16,606	65,431
	2	0.40	0.35	53	90	5	3,726	5,022			0.60	0.55	57	90	5	4,487	7,015
	3	0.60	0.55	55	150	14	5,360	8,379			0.60	0.55	55	150	14	8,662	13,541
	4	0.65	0.60	51	240	17	7,473	11,973			0.55	0.50	54	210	4	3,963	6,021
	5	0.50	0.45	53	330	7	4,338	6,374			0.60	0.55	57	270	13	8,335	13,031
	6	Unused									0.90	0.10	54	330	7	5,814	9,216
2,500	1	0.90	0.55	53	0	23	10,724	18,172	45,287		0.90	0.10	55	30	16	11,215	17,779	70,865
	2	0.65	0.60	52	120	19	8,942	14,328			0.60	0.55	57	90	5	5,002	7,819
	3	0.65	0.60	51	240	17	7,980	12,786			0.60	0.55	55	150	14	9,176	14,346
	4	Unused									0.55	0.50	54	210	4	4,482	6,810
	5	Unused									0.60	0.55	57	270	13	8,850	13,835
	6	Unused									0.65	0.25	54	330	7	6,790	10,276

Table 6. Table 6 : Summary of optimal GPS & Galileo servicing depots from OFLP instances with m s , dry = 1 , 000 subscript 𝑚 𝑠 dry 1 000 m_{s,\text{dry}}=1,000   kg times absent kg \text{\,}\mathrm{k}\mathrm{g}

		$D=1$								$D=2$
$m_{d,\text{dry}}$, $\mathrm{kg}$	No.	$a$, $\mathrm{DU}$	$e$	$i$, $\mathrm{deg}$	$\Omega$, $\mathrm{deg}$	${𝒦}_a$	$m_{d,\text{wet}}$, $\mathrm{kg}$	Depot EMLEO, $\mathrm{kg}$	Total EMLEO, $\mathrm{kg}$		$a$ $\mathrm{DU}$	$e$	$i$ $\mathrm{deg}$	$\Omega$ $\mathrm{deg}$	${𝒦}_a$	$m_{d,\text{wet}}$ $\mathrm{kg}$	Depot EMLEO, $\mathrm{kg}$	Total EMLEO, $\mathrm{kg}$
1,500	1	1.00	0.10	55	30	16	7,028	11,392	47,019		1.00	0.10	55	30	16	11,811	19,146	74,732
	2	0.60	0.55	58	90	5	3,377	5,279			0.60	0.55	58	90	5	5,211	8,146
	3	0.95	0.10	55	150	14	6,208	9,956			0.95	0.10	55	150	14	10,183	16,331
	4	0.60	0.55	53	210	4	2,886	4,513			0.60	0.55	53	210	4	4,230	6,612
	5	0.95	0.15	55	270	13	5936	9,609			1.00	0.10	55	270	13	9,659	15,658
	6	0.95	0.05	55	330	7	3,948	6,270			1.00	0.00	56	330	7	5,564	8,839
2,000	1	0.95	0.10	55	30	16	7,849	12,589	53,192		1.00	0.10	55	30	16	12,559	20,358	81,195
	2	0.60	0.55	58	90	5	3,891	6,084			0.60	0.55	58	90	5	5,725	8,950
	3	0.60	0.55	55	150	14	7,061	11,038			0.95	0.10	55	150	14	10,927	17,525
	4	0.55	0.50	54	210	4	3,498	5,315			0.60	0.55	53	210	4	4,744	7,416
	5	0.60	0.55	56	270	13	6,833	10,682			0.95	0.10	55	270	13	10,507	16,852
	6	0.95	0.05	55	330	7	4,713	7,484			1.00	0.00	56	330	7	6,353	10,093
2,500	1	0.95	0.10	55	30	16	8,594	13,783	58,786		1.00	0.10	55	30	15	12,653	20,510	90,427
	2	0.60	0.55	58	90	5	4,406	6,888			0.60	0.55	58	90	5	6,240	9,754
	3	0.60	0.55	55	150	14	7,575	11,842			0.95	0.10	55	150	14	11,671	18,719
	4	0.55	0.50	54	210	4	4,018	6,104			0.60	0.55	53	210	4	5,258	8,221
	5	0.60	0.55	56	270	13	7,348	11,487			0.95	0.10	55	270	13	11,252	18,046
	6	0.90	0.10	55	330	7	5,477	8,682			1.00	0.05	56	330	8	9,453	15,176

Table 7. Table 7 : Summary of refined GPS & Galileo servicing depots with m s , dry = 500 subscript 𝑚 𝑠 dry 500 m_{s,\text{dry}}=500   kg times absent kg \text{\,}\mathrm{k}\mathrm{g} and D = 1 𝐷 1 D=1

$m_{d,\text{dry}}$ $\mathrm{kg}$	No.	$a$, $\mathrm{DU}$	$e$	$i$, $\mathrm{deg}$	$\Omega$, $\mathrm{deg}$	$m_{d,\text{wet}}$, $\mathrm{kg}$ (% change)	Depot EMLEO, $\mathrm{kg}$ (% change)	Total EMLEO, $\mathrm{kg}$ (% change)
1,500	1	1.0488	0.0325	56.98	19.39	5,037 (-19.41%)	8,136 (-14.09%)	35,621 (-6.90%)
	2	0.5606	0.5375	54.62	90.59	2,728 (-1.10%)	4,187 (-2.89%)
	3	0.5812	0.5661	54.67	139.74	4,571 (-6.01%)	7,106 (-6.19%)
	4	0.5110	0.4936	54.56	202.21	2,461 (-1.62%)	3,663 (-3.59%)
	5	0.5571	0.5009	54.91	260.28	4,579 (-2.24%)	6,976 (-4.69%)
	6	0.6122	0.3192	54.72	324.97	3,676 (-0.31%)	5,552 (-0.90%)
2,000	1	0.9782	0.0391	56.42	19.54	5,883 (-14.45%)	9,376 (-11.09%)	40,397 (-4.70%)
	2	0.5382	0.5187	54.79	90.62	3,278 (-13.67%)	4,966 (-1.10%)
	3	0.5813	0.5661	54.77	139.79	5,074 (-5.64%)	7,889 (-5.85%)
	4	0.6218	0.5829	51.84	243.28	7,501 (+0.38%)	11,871 (-0.86%)
	5	0.4831	0.4770	54.18	325.73	4,306 (-0.74%)	6,295 (-1.24%)
2,500	1	0.8762	0.5199	52.68	1.99	10,781 (+0.53%)	18,116 (-0.31%)	45,008 (-0.62%)
	2	0.6285	0.5871	52.42	120.60	8,967 (+0.27%)	14,232 (-0.67%)
	3	0.6204	0.5818	51.51	243.56	8,004 (+0.30%)	12,659 (-0.99%)

Table 8. Table 8 : Keplerian Elements of GPS Constellation

Satellite number	$a$, $\mathrm{km}$	$e$	$i$, $\mathrm{deg}$	$\Omega$, $\mathrm{deg}$	$\omega$, $\mathrm{deg}$
1	26560.355	6.4584e-03	55.53	150.07	53.20
2	26560.460	4.7800e-03	54.18	72.93	188.43
3	26561.192	1.3721e-02	55.12	146.99	254.56
4	26561.008	1.2823e-02	55.42	267.35	41.45
5	26560.439	2.4678e-02	55.07	17.50	309.60
6	26560.919	8.8526e-03	55.91	328.36	127.48
7	26572.909	2.0378e-02	55.39	17.68	280.51
8	26560.094	1.4080e-02	55.97	325.81	276.13
9	26559.723	1.0584e-02	54.70	203.57	25.15
10	26560.771	8.3765e-03	55.43	266.30	75.05
11	26560.023	1.4504e-02	53.36	134.59	65.84
12	26559.797	2.0014e-03	56.10	326.58	147.45
13	26559.858	1.6622e-02	54.46	202.48	232.26
14	26559.181	5.7239e-03	55.19	79.74	64.40
15	26559.538	1.0562e-02	54.72	261.70	56.85
16	26560.119	1.1835e-02	56.66	23.12	53.36
17	26560.353	1.3191e-02	53.52	197.47	48.46
18	26560.199	1.0589e-02	55.60	322.15	38.13
19	26560.529	6.0448e-03	53.61	203.00	209.29
20	26559.042	2.7915e-03	56.62	22.64	304.45
21	26560.934	2.4195e-03	54.72	141.08	112.11
22	26561.428	4.2503e-03	55.93	82.28	58.76
23	26559.493	7.3252e-03	53.62	258.80	20.82
24	26559.720	7.7127e-03	55.09	320.89	7.69
25	26560.447	8.2834e-03	55.92	82.17	217.98
26	26560.306	6.3642e-03	54.94	141.77	229.63
27	26561.026	2.1567e-03	55.12	144.23	187.36
28	26560.756	2.9509e-03	55.74	23.41	183.78
29	26559.913	2.8304e-03	55.64	80.65	180.09
30	26560.207	2.4235e-03	54.48	264.26	191.16
31	26560.209	8.7880e-04	55.25	25.27	207.44

Table 9. Table 9 : Keplerian Elements of Galileo Constellation

Satellite number	$a$, $\mathrm{km}$	$e$	$i$, $\mathrm{deg}$	$\Omega$, $\mathrm{deg}$	$\omega$, $\mathrm{deg}$
1	29600.198	4.8800e-05	57.04	17.43	2.09
2	29600.168	2.7400e-04	57.04	17.46	298.10
3	29600.354	1.8980e-04	55.18	137.75	231.25
4	29600.327	1.5750e-04	55.18	137.80	138.07
5	27977.498	1.6173e-01	50.18	326.56	129.27
6	27977.430	1.6145e-01	50.21	325.63	130.07
7	29600.058	4.4700e-04	56.84	17.45	239.84
8	29600.127	4.4590e-04	56.84	17.53	226.09
9	29600.445	5.0140e-04	55.54	258.02	2.52
10	29600.444	3.5460e-04	55.54	258.04	345.16
11	29600.271	1.8730e-04	55.20	137.45	344.51
12	29600.294	1.3840e-04	55.19	137.52	1.47
13	29600.440	4.7190e-04	55.68	257.99	324.87
14	29600.265	4.0670e-04	55.68	257.97	296.90
15	29600.433	1.5940e-04	54.86	137.62	310.57
16	29600.438	1.6710e-04	54.85	137.64	29.68
17	29600.412	1.6100e-04	54.86	137.60	248.61
18	29600.434	7.4200e-05	54.85	137.64	230.50
19	29600.355	3.4110e-04	55.66	257.78	293.86
20	29600.347	5.1140e-04	55.66	257.84	290.54
21	29600.349	4.7900e-04	55.66	257.82	302.25
22	29600.344	5.4600e-04	55.66	257.79	285.85
23	29600.073	4.5590e-04	57.22	17.33	259.06
24	29600.076	3.9160e-04	57.22	17.38	263.93
25	29600.077	3.5470e-04	57.22	17.36	262.57
26	29600.077	4.2110e-04	57.22	17.30	260.34
27	29600.059	3.8470e-04	57.21	17.22	213.15
28	29600.059	2.7530e-04	57.21	17.21	204.09

Equations

$$\dot{\mathbf{x}}=\mathbf{B}(\mathbf{x})\mathbf{F}+\mathbf{D}(\mathbf{x})$$

$$\mathbf{F}=[F_r,F_{\theta },F_n{]}^{\mathrm{T}}$$

$$p=a(1-e^2)$$

$$f=e\cos (\Omega +\omega )$$

$$g=e\sin (\Omega +\omega )$$

$$h=\tan \left(\frac{i}{2}\right)\cos \Omega$$

$$k=\tan \left(\frac{i}{2}\right)\sin \Omega$$

$$L=\Omega +\omega +\theta$$

$$\boldsymbol{B}(\boldsymbol{x})=\left[\begin{array}[]{ccc}\dfrac{2a^{2}}{h}e\sin\theta&\dfrac{2a^{2}p}{hr}&0\\[1.00006pt] \sqrt{\dfrac{p}{\mu}}\sin L&\sqrt{\dfrac{p}{\mu}}\dfrac{1}{w}[(w+1)\cos L+f]&-\sqrt{\dfrac{p}{\mu}}\dfrac{g}{w}[h\sin L-k\cos L]\\[1.00006pt] -\sqrt{\dfrac{p}{\mu}}\cos L&\sqrt{\dfrac{p}{\mu}}[(w+1)\sin L+g]&\sqrt{\dfrac{p}{\mu}}\dfrac{f}{w}[h\sin L-k\cos L]\\[1.00006pt] 0&0&\sqrt{\dfrac{p}{\mu}}\dfrac{s^{2}}{2w}\cos L\\[1.00006pt] 0&0&\sqrt{\dfrac{p}{\mu}}\dfrac{s^{2}}{2w}\sin L\\[1.00006pt] 0&0&\sqrt{\dfrac{p}{\mu}}\dfrac{1}{w}[h\sin L-k\cos L]\end{array}\right]$$

$$\boldsymbol{D}(\boldsymbol{x})=\left[\begin{array}[]{llllll}0&0&0&0&0&\sqrt{\mu p}\left(\dfrac{w}{p}\right)^{2}\end{array}\right]^{\mathrm{T}}$$

$$\mathbf{F}=\frac{T}{m}{\left[\cos \alpha \cos \beta ,\sin \alpha \cos \beta ,\sin \beta \right]}^{\mathrm{T}}$$

$$Q=\left(1+W_{p}P\right)\sum _{\text{œ}}{S}_{\text{œ}}{W}_{\text{œ}}{\left(\frac{\text{œ}-{\text{œ}}_T}{{\dot{\text{œ}}}_{xx}}\right)}^2,\text{œ}=a,f,g,h,k$$

$$S_{\text{œ}}=\{\begin{array}{ll}{\left[1+{\left(\frac{\left|a-a_T\right|}{\sigma a_T}\right)}^{\nu }\right]}^{1/\zeta },& \text{œ}=a\\ 1,& \text{ otherwise }\end{array}$$

$$P=\exp \left[k_{r_p}\left(1-\frac{r_p}{r_{p\min }}\right)\right]$$

$$r_p=a(1-e)$$

$$\underset{\alpha ,\beta }{\min }\dot{Q}$$

$$\dot{Q}$$

$$D_1$$

$$D_2$$

$$D_3$$

$${\alpha }^{\ast }$$

$${\beta }^{\ast }$$

$$Z=Z_d{Z}_l=\exp \left(\frac{\Delta V_2}{g_0{I}_{\text{sp},d}}\right)\exp \left(\frac{\Delta V_1}{g_0{I}_{\text{sp},l}}\right)$$

$$r_p=a(1-e),r_a=a(1+e)$$

$$\{\begin{array}{ll}\begin{array}{rl}\Delta V_{p,1}& =\sqrt{\mu \left(\frac{2}{r_0}-\frac{2}{r_0+r_p}\right)}-\sqrt{\frac{\mu }{r_0}}\\ \Delta V_{p,2}& =\sqrt{\mu \left(\frac{2}{r_p}-\frac{1}{a}\right)}-\sqrt{\mu \left(\frac{2}{r_p}-\frac{2}{r_0+r_p}\right)}\end{array}& \text{ if second burn at }{r}_p\\ \begin{array}{rl}\Delta V_{a,1}& =\sqrt{\mu \left(\frac{2}{r_0}-\frac{2}{r_0+r_a}\right)}-\sqrt{\frac{\mu }{r_0}}\\ \Delta V_{a,2}& =\sqrt{\mu \left(\frac{2}{r_a}-\frac{1}{a}\right)}-\sqrt{\mu \left(\frac{2}{r_a}-\frac{2}{r_0+r_a}\right)}\end{array}& \text{ if second burn at }{r}_a\end{array}$$

$$Z=\min \left\{Z_p,Z_a\right\}$$

$$Z_p=(Z_d{Z}_l{)}_p$$

$$Z_a=(Z_d{Z}_l{)}_a$$

$$\underset{X,Y}{\min }$$

$$\sum _{j=1}^n{X}_{ij}=1\forall i\in \mathcal{K}$$

$$\sum _{i=1}^k{X}_{ij}\le Y_j\forall j\in \mathcal{N}$$

$$X_{ij},Y_j\in \{0,1\}\forall i\in \mathcal{K},j\in \mathcal{N}$$

$$\underset{X,Y}{\min }$$

$$m_{d,\text{dry}}{Z}_{d,j}{Y}_j+\sum _{i=1}^k{D}_i({\tilde{\mathbf{c}}}_{ij}+m_{s,L})Z_{d,j}{X}_{ij}\le m_{l,\max }\forall j\in \mathcal{N}$$

$$\text{eqref}eq:FL{P}_{p}ro{b}_{c}onstrain{t}_{o}n{e}_{a}llocation-\text{eqref}eq:FL{P}_{p}ro{b}_{c}onstrain{t}_{b}inary$$

$${\tilde{\mathbf{c}}}_{ij}=m_{s,1}-m_{s,\text{dry}}-m_{s,L}$$

$$m_{s,1}=m_{s,2^{-}}+\Delta m_a=m_{s,2^{+}}+m_{s,L}+\Delta m_a$$

Full text

Nomenclature

1 Introduction

On-orbit servicing, assembly, and manufacturing (OSAM) is a key piece of technology gaining attention from both commercial and governmental players alike [1, 2, 3]. Examples include the Robotic Servicing of Geosynchronous Satellites (RSGS) program by the Defense Advanced Research Projects Agency (DARPA) [4] and Lockeed Martin’s subsidiary SpaceLogistics [5] for servicing geostationary orbit (GEO) satellites, or NASA’s OSAM-1 [6] that focuses on polar low-Earth orbit (LEO). The primary motivation for OSAM has been the cost reduction that comes with extending the lifespan of satellites [7, 8, 9, 10]. It is also possible to recognize the possibility for OSAM to contribute to a more sustainable practice in the space sector, both in terms of the orbital environment and the environmental impact of rocket launch [11].

To date, multiple studies on OSAM applications for GEO and low-inclination geosynchronous orbit (GSO) satellites have been conducted, driven by the relatively high cost of assets lying in these orbits [12, 13, 14, 15]. A particular convenience of GEO servicing comes from the fact that all potential client satellites lie on the equatorial plane; hence, the servicer simply needs to conduct a phasing maneuver to deliver the service. Sarton du Jonchay and Ho [16] studied servicing architectures including both MEO and GEO satellites, but assumed these to be coplanar. Meanwhile, studies for LEO has also been considered by Luu and Hastings [17, 18, 19] and Sirieys et al [20]. Luu and Hastings [18] also provides a review of ongoing work on OSAM for LEO constellations. Converse to GEO constellations, LEO constellations exist across multiple orbital planes, but strong $J_{2}$ perturbations may be leveraged to use differential drift to the servicer’s advantage.

A substantially different scenario from an astrodynamics perspective arises when servicers must navigate between orbital planes of the constellation by primarily utilizing their own propulsion system. Such maneuvers are commonly prohibitively expensive for chemical thrusters, thus requiring the use of low-thrust propulsion with high specific impulse. This scenario is relevant when the constellation to be serviced is at high altitudes, such as Medium Earth Orbit (MEO) or high-inclination GSO. Fortunately, if the servicing is limited to a single constellation, certain features such as the semimajor axis, inclination, eccentricity, and argument of perigee are commonly shared among the constellation fleet. Furthermore, OSAM needs such as refueling and orbit alteration, which are the primary OSAM needs that have recently occurred from both private and governmental players [21, 22, 23, 24], can be planned well in advance. This gives the servicing vehicle enough time to carry out an economical transfer with a longer duration to reach its clients. Hall and Papadopoulous [25] have previously reported on the hardware considerations, while Leisman et al [26] conducted a systems engineering study for servicing the GPS constellation.

The placement of facilities is a fundamental problem in terrestrial logistics. The facility location problem (FLP), since its initial introduction by Cooper [27], has been studied extensively over the past decades [28, 29, 30, 31, 32]. There exist many flavors of FLPs, notably with the distinction between discrete and continuous formulations. In the discrete case, a discrete set of potential locations are considered, with a flexible number of facilities. In contrast, in the continuous case, facility locations are considered in continuous space, but the number of facilities must typically be prescribed [33].

The FLP has seen a wide range of applications, including emergency services [34], telecommunications [35], healthcare [36], waste disposal management [37], and disaster response [38]. In the context of space applications, McKendree and Hall [39] studied a single-facility variant of the FLP to locate an interplanetary manufacturing plant. Dorrington and Olsen [40] applied a variant of the FLP to the asteroid mining problem. Finally, Zhu et al [41] apply the FLP for OSAM depots in Sun-synchronous orbits (SSO).

This work proposes a method to identify optimal combinations of facility locations where the servicing depots are to be placed for an OSAM architecture of a particular satellite constellation. This is done by adopting the FLP to the on-orbit servicing depot of high-altitude satellite constellations; the proposed formulation is coined as the Orbital Facility Location Problem (OFLP). The OFLP is able to simultaneously find the optimal number of facilities, their orbital position within a discretized set of candidate slots, as well as the optimal allocation of each client to the appropriate facility.

The two key modifications of the OFLP come from (1) the formulation of the allocation cost, and (2) the formulation of a facility’s usage cost. As is the case in terrestrial logistics, the optimal location(s) to place facilities depend on some form of a “distance” metric between the facility and the client(s) it must service. Furthermore, the use of an additional facility has associated incurred costs, which must also be taken into account. However, in the context of space-based facilities, the “distance” cannot be simple Euclidean norms or Manhattan distance, but must rather account for the cost of the trajectory to be taken by the servicer. In this work, a Lyapunov feedback controller known as Q-Law [42, 43, 44] is utilized to compute the cost associated with conducting a return trip between a facility and its client satellite using low-thrust propulsion. The total propellant required for these two maneuvers is used as the “distance”. To the best of the authors’ knowledge, this work is the first attempt at considering low-thrust transfers in a FLP framework.

The costs incurred in building a facility at a given location in space are considered on the basis of the launch and orbit insertion cost of the depot carrying with it its servicers, propellant, and payload. In order to coherently consider the orbit transfer costs between the facilities and clients as well as the facilities’ building cost through a single objective function, this work formulates the optimization problem in terms of effective mass to LEO (EMLEO).

The proposed OFLP framework enables the non-trivial, optimal placement of depots for OSAM of spatially distributed clients such as for satellite constellations. By leveraging the discrete version of the FLP, the number of facilities does not need to be predetermined, thus allowing parametric studies in terms of high-level parameters of the architecture, such as the dry mass of the depot or the number of trips to be made to each client.

Once an architecture design is obtained by the OFLP, this work further proposes a refinement scheme whereby the location of each facility is re-optimized in continuous space, while freezing the allocations determined by OFLP. This increases the fidelity of the designed OSAM architecture as the facility location is no longer restricted by the coarseness of the discretized orbital slots.

This paper is organized as follows. In Section 2, the design of low-thrust transfers and the associated approach for estimating the rendez-vous cost for a given pair of orbital elements is introduced. Then, Section 3 discusses the FLP along with its adaptation to the OSAM case. Section 4 presents numerical results applying the formulation to service a combination of the GPS and Galileo constellations. Finally, Section 5 provides a conclusion to this work.

2 Astrodynamics Background

Typical on-orbit servicing demands are well-known in advance; this results in phasing maneuvers having little impact on overall transfer costs, as arbitrarily long phasing maneuver time may be utilized to minimize propellant expenditure. As such, when considering the transfer cost in the on-orbit servicing context, the primary propellant expenditure comes from the orbit transfer cost, i.e. the cost associated with modifying the orientation and shape of the orbit in 3D space. To this end, a Lyapunov feedback controller for spacecraft trajectories, commonly referred to in the astrodynamics literature as Q-law, is employed.

2.1 Orbit Transfer via Q-Law

The orbital transfer is conducted using Q-law, a Lyapunov-function-based feedback control law first introduced by Petropoulos [42, 43, 44]. Despite its sub-optimal nature, its ability to rapidly compute multi-revolution low-thrust transfers has made it a popular tool, particularly for large-scale preliminary transfer designs [45, 46, 47, 48, 49]. The feedback law consists of determining the optimal thrust direction given the current and target orbital elements. Among the target elements, only the five slow variables are used for the transfer, as the considered problem does not necessitate rendez-vous.

The dynamics of the satellite may be expressed in terms of orbital elements via Gauss’ variational equations, or the Variation of Parameters (VOP) equations. Let the state be given by $\boldsymbol{x}\in\mathbb{R}^{6}$, and consider $\boldsymbol{B}\in\mathbb{R}^{6\times 3}$ and $\boldsymbol{D}\in\mathbb{R}^{6}$, such that

[TABLE]

where $\boldsymbol{F}\in\mathbb{R}^{3}$ is the perturbing force expressed in terms of radial, tangential, and normal components

[TABLE]

The set of elements $\boldsymbol{x}$ may be Keplerian elements $[a,e,i,\Omega,\omega,\theta]$ or other alternative sets of elements. While the original Q-law has been developed in terms of Keplerian elements, the well-known singularities at $i=0$ and $e=0$ as these are typical orbital elements in which a spacecraft may be placed. For this reason, the use of alternative element sets, such as the modified equinoctial elements (MEE), given in terms of Keplerian elements as

[TABLE]

is beneficial as the singularity is moved to $i=\pi$. Furthermore, the use of the MEE with the semi-parameter $p$ replaced by the semimajor axis $a$ has been previously reported to yield convergence benefits [50], and is employed in this work as well. For this set of elements $\boldsymbol{x}=[a,f,g,h,k,L]$, the VOP are given by

[TABLE]

and

[TABLE]

Note that $\boldsymbol{F}$ can be due to any form of perturbing acceleration, such as propulsive force, atmospheric drag, third-body attraction, or non-spherical gravity. The acceleration due to the propulsive force, which is to be determined to guide the spacecraft to its target orbit, is given by

[TABLE]

where $\alpha$ and $\beta$ are the in-plane and out-of-plane angles, respectively.

Denoting the osculating element $\text{\oe}\in[a,f,g,h,k]$ and the targeted elements as $\text{\oe}_{T}$, the Lyapunov function is given by

[TABLE]

In essence, $Q$ penalizes the difference between œ and $\text{\oe}_{t}$ through the subtraction term in the summation. $S_{\text{\oe}}$ is a scaling factor given by

[TABLE]

where $\sigma$, $\nu$ and $\zeta$ are scalar coefficients, which prevents non-convergence of $a\to\infty$. This is necessary as when $a\to\infty$, the $\dot{\text{\oe}}_{xx}$ terms also tend to $\infty$ and thus $Q$ is reduced, however, this is not a physically useful solution. $W_{\text{\oe}}$ are scalar weights that may be assigned to different elements, if one is to be favored for targeting over another. $P$ is a penalty term on the periapsis radius, given by

[TABLE]

where $r_{p}$ is the current orbit’s periapsis radius given by

[TABLE]

and $r_{p\min}$ is a user-defined threshold. Here, $k_{r_{p}}$ is also a pre-defined constant on this penalty term and represents the gradient of this exponential barrier function near $r_{p\min}$. $W_{p}$ is a scalar weight to be placed on the periapsis penalty term. The $\dot{\text{\oe}}_{xx}$ terms represent the maximum rates of change of a given orbital element with respect to the thrust direction and true anomaly along the osculating orbit and are given in the Appendix.

Through the application of Lyapunov control theory, the Q-law strategy consists of choosing the control angles $\alpha$ and $\beta$ such that the time-rate of change of $Q$ is minimized at each time-step

[TABLE]

where $\dot{Q}$ can be expressed using the chain rule as

[TABLE]

where

[TABLE]

The choice of $\alpha$ and $\beta$ based on condition (11), given by

[TABLE]

ensures the fastest possible decrease of $Q$, thereby providing the best action for the spacecraft to take to arrive at $\text{\oe}_{t}$. Note that while $\dot{\text{\oe}}$ consists simply of the first 5 rows of the VOP given in expression (4), the expression for $\frac{\partial Q}{\partial\text{\oe}}$ is cumbersome to derive analytically. Instead, a symbolic toolbox is used to obtain these expressions. Finally, it is noted that there are several work in the past that extended the Q-law scheme to include a mechanism for coasting [44, 48] as well as other perturbing forces such as oblateness or three-body effects [50, 45, 47]. While these variants are not considered in this work, they may be used in place of the basic controller in the OFLP, which will be presented subsequently.

2.2 Launch and Insertion Costs of Depot

The launch and insertion cost of a given depot depends on the mass of the depot and its servicer(s) as well as its orbital elements. In the context of constellation servicing, the depots are to be placed in an orbit sharing the same inclination as the constellation fleet. Hence, out of the six Keplerian elements, only the semimajor axis and eccentricity have an effect on the launch and insertion cost.

In order to place the depot into its desired orbital slot, the launch vehicle must first lift off the depot from the launch pad and place it into a transfer orbit, and the depot must then conduct an insertion maneuver upon arrival at its desired orbit. Both the insertion cost into the transfer orbit and the insertion cost to the depot’s desired orbit depend on the choice of the orbital slot; as such, this work considers the effective mass to LEO (EMLEO) as the standardized metric for evaluating the depot’s equivalent, total insertion cost as a standardized mass. EMLEO is the hypothetical mass the depot would have if it were to first be placed into a circular LEO on the same orbital plane as the transfer orbit, rather than into the transfer orbit itself, by the launch vehicle. It thus gives a consistent cost to compare depots that have differing transfer orbit insertion costs.

In this work, we assume a coplanar Hohmann transfer to deliver a depot. With the aforementioned definition, the EMLEO of a given orbital slot can be found by calculating the propellant mass needed to perform a Hohmann transfer from a (hypothetical) LEO with radius $r_{0}$ to the desired orbital slot. The first impulse of the Hohmann transfer corresponds to the transfer orbit insertion, which is (hypothetically) carried out by the launch vehicle, while the second impulse of the Hohmann transfer is carried out by the depot itself. Using the rocket equation, the combined mass-ratio $Z$ is the product of the mass-ratio due to the first burn by the launch vehicle $Z_{l}$ and the mass-ratio due to the second burn by the depot itself $Z_{d}$, and can be expressed as

[TABLE]

where $\Delta V_{1}$ and $\Delta V_{2}$ are the impulsive burn magnitudes of the two maneuvers, $g_{0}$ is the standard acceleration due to gravity, $I_{\text{sp},d}$ is the specific impulse of the depot, and $I_{\text{sp},l}$ is the launch vehicle. To compute the $\Delta V$ magnitudes, the perigee, and apogee are first computed from the Keplerian elements of the $j^{\mathrm{th}}$ depot via

[TABLE]

To transfer from a circular orbit to an elliptical orbit via a Hohmann transfer, the elliptical transfer orbit would have a periapsis at the radius of the initial circular orbit and an apoapsis at either the final orbit’s perigee or apogee. Since $I_{\text{sp},d}$ and $I_{\text{sp},l}$ are not necessarily the same, both options are considered, and the strategy resulting in a smaller mass-ratio $Z_{j}$ is chosen. In each case, the $\Delta V$’s are given by

[TABLE]

Then, the mass-ratio $Z$ is the combination of $Z_{d}$ and $Z_{l}$ that results in a smaller value

[TABLE]

where

[TABLE]

Finally, note that since both $Z_{d}$ and $Z_{l}$ are ratios that do not depend on any mass; hence, $Z$ can be pre-computed once for a given orbital slot and reused regardless of the facility’s mass. By multiplying the depot’s mass at the start of its operation after the arrival insertion burn with $Z$, the EMLEO of the depot is obtained.

3 Methods

This section first introduces the static facility location problem (FLP), and discusses its adoption for the orbital facility case, resulting in the orbital facility location problem (OFLP). In doing so, discretized candidate locations in space for the facilities must be defined. Considerations involved in this process are discussed in this section. Then, the cost coefficients of the FLP must be formulated for the orbital facility case; specifically, the allocation cost consisting of the required propellant for a round-trip between a facility and its client is introduced. Finally, once the OFLP is solved, the orbital location of each facility is refined in continuous space, while maintaining the allocations of clients determined by the OFLP. Figure 1 illustrates the flow of the proposed method.

3.1 Overview of Facility Location Problem

Given a set of clients $\mathcal{K}=\{1,2,\ldots,k\}$ and a set of potential sites for facilities $\mathcal{N}=\{1,2,\ldots,n\}$, the single-source capacitated FLP aims to locate $k\leq n$ facilities to service all $m$ clients in a cost-optimal sense. The problem is a binary linear program (BLP), given by

[TABLE]

The decision variables $X\in\mathbb{B}^{k\times n}$ and $Y\in\mathbb{B}^{n}$ are restricted to be binary by (20d). The variable $Y_{j}$ dictates whether the facility $j$ is used within the architecture, and the variable $X_{ij}$ dictates whether facility $j$ is assigned to client $i$. The objective of this problem, as given in (20a), is to minimize the sum of costs associated with the allocation of the facility $j$ to client $i$, and the costs of using facility $j$. The entries in the matrix $\boldsymbol{c}_{ij}$ are weights associated with the $i$-$j$ allocation, while the entries of the vector $\boldsymbol{f}_{j}$ are weights associated with using facility $j$. Constraints (20b) ensure that exactly one facility is assigned to each client. Finally, constraints (20c) ensure that the allocation(s) of facility $j$ to client(s) do(es) not exceed the capacity of the facility.

A particular advantage of this formulation of the FLP, compared to its continuous analogs, is that the number of facilities does not need to be known a priori. Rather, the number of facilities is simply given by $\sum_{j=1}^{k}Y_{j}$.

3.2 Orbital Facility Location Problem

The FLP is adapted to the case of on-orbit servicing for satellite constellations, called the Orbital Facility Location Problem (OFLP). The OFLP is formulated as follows:

[TABLE]

Here, $m_{d,\text{dry}}$ is the dry mass of the depot, $m_{s,L}$ is the payload mass to be delivered upon each trip to the client satellites, and $D_{i}$ is the number of trips to be made to each client over the operation period of the depot. $Z_{j}$ is the mass-ratio of launching the depot to the facility site $j$, computed by equation (18), and the entries of $\tilde{\boldsymbol{c}}_{ij}$ are the propellant mass required to conduct a return trip from facility $j$ to deliver the payload to client $i$. The resulting objective (21a) corresponds to the sum of the EMLEO of each facility.

The OFLP takes the same constraints (20b) - (20d) as the original FLP, with an additional constraint (21b) that ensures the wet mass of a given depot does not exceed the maximum launch mass of the launch vehicle $m_{l,\max}$. Note that the wet mass of the depot is computed by multiplying only the mass ratio of the depot’s burn, $Z_{d}$, while the EMLEO in the objective is obtained by multiplying both mass ratios $Z=Z_{d}Z_{l}$.

In effect, the OFLP formulation optimizes both the number of facilities as well as their spatial configuration simultaneously through a cost-metric that is standardized in terms of mass. The choice of the spatial configuration of a facility results from the trade-off between the cost of establishing it into a particular orbit and the total costs of accessing the allocated clients from it.

3.3 Discretization of Facility Location Candidates

In the FLP, candidate facility locations must be represented as a discrete number of options in set $N$. As such, when formulating the OFLP, it is necessary to discretize locations along an orbit into slots. To avoid combinatorial explosion, it is necessary to reduce the dimensions of the discretized space. Firstly, the anomaly-like variable of the facility is omitted; as discussed earlier, this is a common assumption for multi-revolution low-thrust transfers (see more details in the Q-law literature [42, 43, 44]). In our case, the transfers between facilities and clients involve changing any combination of the energy, orbital shape, and orbital plane, within which phasing maneuvers are comparatively inexpensive, especially if the phasing maneuvers can take place over a long time horizon (e.g., when refueling demands for a client satellite can be predicted months in advance). As such, the transfer cost between the facilities and the clients may be modeled as an orbital transfer rather than a rendez-vous problem. Second, if the constellation shares the same inclination among all of its fleet members, the facility’s inclination may be set to match this value. Finally, Earth-based constellations are typically near-circular. Then, under the assumption of pure two-body dynamics, the transfer cost from a facility orbit to a client’s orbit is negligibly affected by the argument of perigee of the initial departure orbit.

Hence, making appropriate assumptions that depend on the distribution of the clients, the dimension of the discretization would vary between three and five, among the semimajor axis, eccentricity, inclination, RAAN, and argument of perigee. For example in Zhu et al [41], the client satellites are distributed on the same inclination and are near-circular; thus, only the semimajor axis, eccentricity, and RAAN have been considered. In the numerical example studied in this work that will be introduced in Section 4, the client satellites are near-circular but on varying inclinations; therefore, the facility locations are discretized in terms of semimajor axis, eccentricity, inclination, and RAAN.

3.4 Computation of Allocation Cost Matrix

The allocation cost matrix $\tilde{\boldsymbol{c}}_{ij}$ has as entries the propellant mass required to allocate client $i$ to facility $j$. This mass is computed by considering the return trip consisting of (a) the outbound trip traveling from the facility to the client satellite, (b) depositing the payload to the client, and (c) the inbound trip returning from the client satellite back to the facility, as illustrated in Figure 2. Since the overall mass of the spacecraft affects the attainable acceleration and consequently the trajectory as well, the propellant mass expenditures are computed backward in time.

[TABLE]

where $m_{s,L}$ is a constant in the problem, and

[TABLE]

The propellant expenditures in (a) and (c) is computed through Q-law

[TABLE]

where $\Delta t_{\text{OT},{ji}}$ is the time taken by the Q-law controller to transfer from facility $j$ to client $i$, $\Delta t_{\text{OT},{ij}}$ is the time taken by the Q-law controller to transfer from client $i$ to facility $j$. Note that $\Delta t_{\text{OT},ji}\neq\Delta t_{\text{OT},ij}$ since the mass of the spacecraft during phase (a) and (c) are different.

It is also possible that the return trip from a particular facility slot to a client is infeasible due to a time of flight that is prohibitively long. In such case, the entry $\tilde{\boldsymbol{c}}_{ij}$ can simply be set to a high value. Alternatively, it is also possible to eliminate the variable $X_{ij}$ by freezing its value to [math]. This has the added advantage of reducing the model size and is therefore done in this work.

3.5 Refinement of Facility Location in Continuous Space

In OFLP, the orbital slots had to be discretized to efficiently optimize the number of facilities and the allocation of clients to each of them. Thus, given the number of facilities and the allocation of clients obtained from the OFLP solution, we refine the facility locations further in continuous space by formulating a continuous nonlinear programming problem (NLP). Unlike OFLP, the refinement problem is solved for each facility separately. Specifically, by collecting the variables dictating the location of a facility into a vector $\aleph$, the refinement problem aims at adjusting the facility’s location to minimize its EMLEO, and is given by

[TABLE]

where $\mathcal{K}_{a}$ is the set of indices of clients that are allocated to the given facility. Note also that now, the transfer costs from the facility to each client $\tilde{c}_{i}$ as well as the facility’s mass ratios $Z$ and $Z_{d}$ are all functions of variables $\aleph$, and must therefore be computed online. This online computation of $\tilde{c}_{i}$ involves calling Q-law during the optimization, which may result in numerical difficulties if gradient-based algorithms are applied to solve (25). To circumvent this issue, gradient-free metaheuristic algorithms are employed in this work.

4 Numerical Results

The proposed formulation is applied to a scenario focused on simultaneously servicing the GPS and Galileo constellations. Both the GPS and Galileo constellations are located on Medium Earth Orbits (MEO), at similar ranges of inclination around $53^{\circ}$ to $57^{\circ}$, and a radius of about $26\,560\text{\,}\mathrm{k}\mathrm{m}$ and $29\,600\text{\,}\mathrm{k}\mathrm{m}$, respectively. The exact orbital elements of the two fleets are given in the Appendix, in Tables 8 and 9. The GPS constellation consists of 31 satellites on six orbital planes, and the Galileo constellation [51] consists of 28 satellites on four orbital planes, resulting in 59 client satellites in total. Orbital elements of the GPS and Galileo fleets used in this work are given in the Appendix. Figure 3 shows the distribution of the clients in inclination and RAAN.

The OFLP is implemented in Julia and solved with CPLEX 22.1 [52]. Table 2 lists the problem parameters used in the numerical experiments. The launch vehicle’s maximum payload mass capability of $12950\,\mathrm{kg}$ corresponds to the sub-GEO transfer orbit capability of the Ariane 64 [53], used here as an example launch vehicle constraint. The servicer vehicle’s thrust and $I_{\mathrm{sp},s}$ are taken from Sarton du Jonchay et al [14]. Parameters for Q-Law are based on the work by Petropoulos [44].

4.1 Depots Placements via Orbital Facility Location Problem

The orbital slots used for the GPS and Galileo case are discretized in the four slow Keplerian elements excluding the argument of perigee, $[a,e,i,\Omega]$, as summarized in Table 3, resulting in 23,868 slots. The problem has dimensions $X\in\mathbb{B}^{59\times 23,868}$ and $Y\in\mathbb{B}^{23,868}$, resulting in 1,432,080 variables. As previously mentioned, the $\omega$ is not discretized because both the GPS and Galileo satellites are at very low eccentricities. While a constant $\omega=0\,\mathrm{deg}$ is assumed for all orbital slots, the resulting solution would be unaffected for any other choice of $\omega$, assuming the transfer cost does not change either.

Figure 4 shows the contour of facility mass-ratio to $Z$ for the orbital slots in the semimajor axis-eccentricity space, and the contour of allocation cost $\tilde{c}_{ij}$, corresponding to the servicer propellant mass, from the discretized orbital slots to a hypothetical client lying on the same orbital plane, on a circular orbit with a semimajor axis of $1.0$ $\text{\,}\mathrm{D}\mathrm{U}$. The missing entries in $\tilde{c}_{ij}$ are due to locations resulting in infeasible transfers for $r_{p\min}=6878$ $\text{\,}\mathrm{k}\mathrm{m}$. As expected, the orbital slots closer to the facility yield a lower servicer propellant mass, but tend to come at a larger $Z$, which translates to a larger launch cost. The trade-off between these two factors is taken into account by the OFLP.

The effect of varying the demand $D$, the servicer dry mass $m_{s,\text{dry}}$, and the depot dry mass $m_{d,\text{dry}}$ are investigated by constructing and solving OFLP instances for all possible combinations of these parameters. Table 4 shows the performance of CPLEX on a computer with 8 core i7-10700 CPU and 32 GB computer memory. Tables 5 and 6 provide the orbital elements as well as the wet mass and number of clients allocated to each facility. Figures 5 and 6 show the distribution of facilities in semimajor axis-RAAN space. Figures 7 and 8 show the orbits of the facilities in each on their perifocal frame rotated about the $z$-axis by $\Omega$.

Some intuitions can be obtained from these results. Firstly, as expected, facilities tend to align with the clusters of clients, located at 6 distinct regions in $\Omega$-space. The slight misalignment between the facility and the centroid of the cluster is an artifact of the discretization of orbital slots.

Also, as it is visible in both Figures 7 and 8, orbits of facilities can be categorized into ellipses with low eccentricities of around $0.2$ to $0.35$ and ellipses with higher eccentricities of $0.5$ or higher. For the same semimajor axis, the former type requires a larger facility location mass ratio $Z$ but the servicer transfer cost is lower, while the latter type has a lower $Z$ and higher servicer transfer cost. Notably, facilities servicing multiple locations in RAAN, namely facility 4 in $D=1$ and $m_{d,\text{dry}}=2,000$ $\text{\,}\mathrm{k}\mathrm{g}$ as well as all three facilities in $D=1$ and $m_{d,\text{dry}}=2,500$ $\text{\,}\mathrm{k}\mathrm{g}$, have relatively high eccentricities.

Finally, Figure 6 immediately reveals a noteworthy allocation for $D=2$ and $m_{d,\text{dry}}=2,500$ $\text{\,}\mathrm{k}\mathrm{g}$ to facility 6 at $\Omega=330$ $\text{\,}\mathrm{d}\mathrm{e}\mathrm{g}$, stretching a client in the cluster close to $\Omega=30$ $\text{\,}\mathrm{d}\mathrm{e}\mathrm{g}$. This is due to the maximum launch mass constraint (21b) which restricts the wet mass of the facility from exceeding $m_{l,\max}$. Indeed, facility 1 at $\Omega=30$ $\text{\,}\mathrm{d}\mathrm{e}\mathrm{g}$ has a wet mass of $12,653$ $\text{\,}\mathrm{k}\mathrm{g}$, which is very close to $m_{l,\max}=12,950$ $\text{\,}\mathrm{k}\mathrm{g}$. As a result, facility 6 takes up an additional $8^{\mathrm{th}}$ client, as seen in Table 6. It is also noted that this solution with a divided allocation of clients in a cluster to multiple facilities is more challenging to CPLEX, and requires a longer solve time of approximately 45 minutes, compared to the other solutions which are found in up to around 10 minutes, as shown in Table 4.

In the following subsections, the effect of the three parameters, $D$, $m_{s,\text{dry}}$, and $m_{d,\text{dry}}$, is analyzed in further detail.

4.1.1 Effect of Facility Dry Mass

An increase in facility dry mass $m_{d,\text{dry}}$ means an increase in facility launch cost; two things are noted with respect to increasing $m_{d,\text{dry}}$. Firstly, looking at facilities located at the same RAAN, an increase in $m_{d,\text{dry}}$ is seen to result in some cases to a decrease in $a$ and an increase in $e$, provided that the same clients are still allocated to the specific depot. This is due to the lower $Z$ that can be achieved when increasing $e$, as visible from the contour of $Z$ in Figure 4. While decreasing $a$ would also decrease $Z$, decreasing the energy of the facility leads to a more significant penalty on the transfer cost $\tilde{c}_{ij}$ than the penalty associated with increasing $e$ of the facility. This is seen for example in Figure 5 for the depot at $\Omega=30\,\mathrm{deg}$ for $D=1$, or the depot at $\Omega=330\,\mathrm{deg}$ with $D=2$. However, more commonly, the depot locations are observed to be mostly unaffected by an increase in $m_{d,\text{dry}}$. This variability of behavior can be attributed to the discretization of the facility slots; if the benefit of moving a depot is too small to be captured by adjacent facility slots, the optimal solution would not move the depot. It can thus be discerned that the corresponding trend is relatively weak for the chosen coarseness of the facility slots.

Secondly, an increase in $m_{d,\text{dry}}$ is also observed to occasionally alter the distribution of the clients among the depots as well. Optimal solutions bundling two clusters of clients into a single facility begin to appear, as long as the facility mass does not exceed the maximum launch mass, and the required number of return trips dictated by $D$ does not increase. This is observed for $D=1$, where the clusters of clients in the fourth quadrant of $\Omega$ are bundled when increasing $m_{d,\text{dry}}$ from $1.500\text{\,}\mathrm{k}\mathrm{g}$ to $2.000\text{\,}\mathrm{k}\mathrm{g}$, and the other four clusters are also bundled when increasing $m_{d,\text{dry}}$ further to $2.500\text{\,}\mathrm{k}\mathrm{g}$.

It is also noteworthy that in these cases, a facility that services multiple clusters of clients is consistently located at an inclination of $53$ $\text{\,}\mathrm{d}\mathrm{e}\mathrm{g}$ or lower, despite the fact that all clients except for one have inclinations of $53$ $\text{\,}\mathrm{d}\mathrm{e}\mathrm{g}$ or higher. This is a result of the transfer cost involved when the servicer has to make significant changes in RAAN during the round-trip; a lower inclination is more beneficial for such a maneuver providing a larger “lever-arm” to rotate the orbital plane, even though there will be added cost for having to adjust the inclination to match the clients’.

4.1.2 Effect of Demand

The demand $D$ reflects the number of round trips to be conducted between a depot and its clients. While this does not affect the cost of the transfer itself, the additional round trip requires additional propellant and payload to be brought and stored at the depot. As such, as $D$ increases, the saving in propellant mass resulting from placing the depot closer to the client increases. For the depots at the same $\Omega$ between $D=1$ and $D=2$, the orbits of the facilities are generally moved “closer” to the clients, by either increasing $a$, reducing $e$, or both. Examples of such depots include the one at $\Omega=30$ $\text{\,}\mathrm{d}\mathrm{e}\mathrm{g}$ in all cases except for $m_{d,\text{dry}}=1,500$ $\text{\,}\mathrm{k}\mathrm{g}$ and $m_{s,\text{dry}}=1,000$ $\text{\,}\mathrm{k}\mathrm{g}$, as presented in Tables 5 and 6. The existence of some depots remaining in the slots despite the increase in $D$ can again be explained by the coarseness of the discretization of the facility slots.

4.1.3 Effect of Servicer Dry Mass

Increasing the servicer dry mass $m_{s,\text{dry}}$ increases the required propellant mass because the acceleration of the servicer decreases for a given propulsion system of fixed thrust, resulting in a longer transfer. As such, similarly to the effect of $D$, an increase in $m_{s,\text{dry}}$ would move the design towards bringing the depots “closer” to the clients. This effect is most pronounced for facilities servicing a larger number of clients; when the optimal solution uses 6 depots, the facilities located at $\Omega=30$, $150$, and $270$ $\text{\,}\mathrm{d}\mathrm{e}\mathrm{g}$ servicing 16, 14, and 13 clients respectively tend to reside at higher $a$ and lower $e$ comparing the solutions from $m_{s,\text{dry}}=500$ $\text{\,}\mathrm{k}\mathrm{g}$ in Table 5 to the solutions from $m_{s,\text{dry}}=1,000$ $\text{\,}\mathrm{k}\mathrm{g}$ in Table 5.

4.2 Depots Location Refinements

Fixing the number of depots and the allocation of clients to those obtained from the OFLP solution, the depots’ locations are refined in continuous space using metaheuristics. The refinement decision vector consists of the four orbital elements that have been discretized for OFLP, $\aleph=[a,e,i,\Omega]$, and this step is performed for each depot separately. For this demonstration, the facility locations for the three OFLP results with $D=1$ and $m_{s,\text{dry}}=500\,\mathrm{kg}$ are refined. The optimization is done using Differential Evolution [54], with a population size of 50 with a mutation scale factor $F=0.9$.

Table 7 shows the refined locations of the facilities for $m_{s,\text{dry}}=500$ $\text{\,}\mathrm{k}\mathrm{g}$ and $D=1$. The refined results include both facilities that have been refined within the grid of the OFLP and ones that have moved outside in terms of $a$ and $e$. On the other hand, refinements in $i$ have mostly been within the grid except for facility 1 in the $m_{d,\text{dry}}=1,500\,\mathrm{kg}$ and $m_{d,\text{dry}}=2,000\,\mathrm{kg}$ cases, and the refinements in $\Omega$ have strictly been within the grid, to optimally align each facility with the cluster of clients. The occasional fluctuation in $a$, $e$, and $i$ beyond the grid is an artifact of the discretization, where previously in OFLP the transfer costs have been traded-off up to the precision limited by the coarseness of the grid as shown in Figure 4. The resulting reduction in total EMLEO is within 7% for all tested cases. For individual depots, the largest reductions are observed for facility 1 in $m_{d,\text{dry}}=1,500\,\mathrm{kg}$ and $m_{d,\text{dry}}=2,000\,\mathrm{kg}$, where the reduction is within 15%. To understand this fluctuation, taking this facility 1 as an example, Figure 9 shows the contour of the objective (25a) in semimajor axis and eccentricity space, while freezing the inclination and RAAN to the values found by the OFLP or by the refinement NLP, respectively. It is possible to observe that the local minimum of EMLEO clearly shifts as the inclination and RAAN are adjusted, resulting in a noticeably different facility orbit. Overall, the results demonstrate the value of the refinement step by further improving the OFLP results.

Note that these refined solutions are assuming a fixed servicing allocation based on the OFLP results, and so could not have been obtained without first formulating and solving an OFLP instance. Specifically, concurrently optimizing the number of facilities, their locations, and the allocation of each client to each specific facility altogether is prohibitively difficult in a traditional nonlinear global optimization method because it would require large numbers of both binary and continuous variables, which would severely degrade the performance of any global optimization algorithm. In contrast, the proposed scheme combining the OFLP and the continuous refinement problem leverages the strength of both formulations, enabling the design of complex OSAM depot placements.

For the refined depot location solutions, Fig. 10 shows the time of flight and propellant mass of the outbound and inbound transfers. As previously highlighted, the outbound trip requires a longer time and more propellant due to the additional mass to be carried. Also, the time of flight and propellant mass consumption are proportional, as the transfers are assumed to use maximum thrust at all times. Across the three cases, there are roughly three groups of transfers that can be identified:

1. (a)

The first type is the cheapest and fastest, with a time of flight of around 2 days, with propellant consumption of less than 20 $\text{\,}\mathrm{k}\mathrm{g}$; this corresponds to transfers from facilities with a relatively large semimajor axis and low eccentricity servicing clients in close-by orbital planes, such as facility 1 for $D=1$ with $m_{d,\text{dry}}=1,500$ $\text{\,}\mathrm{k}\mathrm{g}$ and $m_{d,\text{dry}}=2,000$ $\text{\,}\mathrm{k}\mathrm{g}$. 2. (b)

The second type of transfer has a time of flight of around 7 to 20 days, and a propellant consumption of around $50$ to $100$ $\text{\,}\mathrm{k}\mathrm{g}$; this corresponds to transfers from facilities with a relatively low semimajor axis and high eccentricity, also servicing clients in similar orbital planes. 3. (c)

The third type of transfer has higher times of flight as well as propellant consumption than the previous two types; this corresponds to transfers to clients in significantly different RAAN.

Figure 11 shows examples of each of these types of transfer for the outbound leg. Type (a) has a time of flight of $1.17\text{\,}\mathrm{d}\mathrm{a}\mathrm{y}\mathrm{s}$ with a propellant consumption of $9.99\text{\,}\mathrm{k}\mathrm{g}$, type (b) has a time of flight of $7.93\text{\,}\mathrm{d}\mathrm{a}\mathrm{y}\mathrm{s}$ with a propellant consumption of $67.95\text{\,}\mathrm{k}\mathrm{g}$, and type (c) has a time of flight of $26.86\text{\,}\mathrm{d}\mathrm{a}\mathrm{y}\mathrm{s}$ with a propellant consumption of $230.06\text{\,}\mathrm{k}\mathrm{g}$.

Noting the GPS Block III launch mass of $3.880\text{\,}\mathrm{k}\mathrm{g}$ or the Galileo launch mass of $700\text{\,}\mathrm{k}\mathrm{g}$ the mass of the depots are of similar scales to a single GPS satellite, and up to a few times heavier than a single Galileo satellite. However, a single facility is able to service and hence extend the lifetime of multiple satellites, hence reducing the overall cost of the programs. Taking for example the servicing architecture with $m_{d,\text{dry}}=1,500\,\mathrm{kg}$ and $D=1$, the total architecture mass is $23,052$ $\text{\,}\mathrm{k}\mathrm{g}$, while replacing all 31 GPS satellites and all 28 Galileo satellites would require launching $139,880$ $\text{\,}\mathrm{k}\mathrm{g}$.

5 Conclusion

This work studies the optimal placement of on-orbit servicing depots for satellite constellations at MEO or higher altitudes. The number of servicing depots as well as their placement(s), from which servicing vehicles fly to rendez-vous with the client satellites, are simultaneously optimized by modifying the facility location problem to formulate the Orbital Facility Location Problem (OFLP). Potential locations for placing the depots are considered in a discretized design space consisting of orbital slots. The combinatorial nature of the orbital slots results in large numbers of transfers to be computed; this is done by using the Q-law backward in time to obtain the return trip propellant mass expenditure. The cost associated with assigning a facility on a particular orbital slot to a client satellite as well as the cost of establishing this facility is combined to form a single objective based on the sum of each facility’s EMLEO. Once a solution for the OFLP has been obtained, the facility locations have been refined using a continuous global optimization scheme while keeping the allocations of clients to each facility fixed to that of the OFLP solution.

The resulting OFLP still retains the form of a binary linear program, which can be solved efficiently even with a large number of variables. The formulation is applied to a scenario for placing depots to service the GPS and Galileo constellations. Multiple instances of the OFLP have been solved with different values for the depot and servicer parameters, and the significant variation in number, as well as the locations of the depots, have been observed. Intuitive design features, such as allocating clients in similar orbital planes to a single depot, have been confirmed to be advantageous by the OFLP, while less intuitive results, such as placing the depot at a slightly lower inclination than the clients to reduce the propellant cost of traversing in RAAN, have also been discovered.

The methods presented in the work, along with the analyses based on the GPS and Galileo satellites, enable studying the feasibility of on-orbit servicing architectures of high-altitude satellite constellations. As on-orbit servicing technology matures in the GEO market, the MEO and GSO markets would present a natural extension of such service, and high-level design of such architecture may be conducted with the OFLP.

Appendix

5.1 Expressions for Q-Law

5.1.1 Maximum rates of change of orbital element $\dot{\text{\oe}}_{xx}$

[TABLE]

where $F$ is the magnitude of the perturbing force.

5.2 Keplerian Elements of Considered Constellations

Tables 8 and 9 show the Keplerian elements of the GPS and Galileo constellation fleet, respectively.