Relative Orbit Control Algorithms and Scenarios for the Inertial Alignment Hold Demonstration Mission by CubeSat Formation Flying

Abstract

CANYVAL-C is a formation-flying mission that demonstrates a coronagraph utilizing two CubeSats. The coronagraph is a space telescope that blocks sunlight to examine the overcast regions around the sun. It is composed of optical and occult segments. Two spacecraft were aligned with respect to an inertial system to configure a virtual telescope using inertial alignment hold technology. The relative orbit control scenario for this mission involves rendezvous, differential air drag control, and inertial alignment holding. Orbit control algorithms and simple strategies that can be automatically constructed onboard have also been developed. For each maneuver, the control performance under the errors from navigation, attitude determination and control, and propulsion systems were assessed via Monte Carlo simulation, taking into account the hardware specifications and operations. In addition to the algorithm and strategy of this mission, the relative orbit control scenario was evaluated for its practicability using Monte Carlo simulations. The feasibility of this mission is ensured by a statistical analysis of the prospect of its success during its operation.

1. Introduction

1.1. Historical Background

Recent advancements in precision flight technology have enabled the operation of numerous formation-flying missions in low-Earth orbit as part of sophisticated distributed satellite systems. Distributed satellite systems have advantages over monolithic architectures in terms of financial efficiency, enhanced flexibility in configuration, and robustness against failure and decommission [1]. The Gravity Recovery and Climate Experiment (GRACE) is aimed at gravimetry by measuring the changes in formation induced by gravitational forces [2]. TSX/TDX (TerraSAR-X/TanDEM-X) is an Earth observation mission that generates a high-precision digital elevation model (DEM) around the globe by achieving a baseline of a few kilometers [3,4]. Prototype Research Instruments and Space Mission Technology Advancement (PRISMA) is a formation-flying technology demonstration mission that investigates model guidance, navigation, and control (GNC) algorithms and hardware [5].

The number of formation-flying missions utilizing nano- and microsatellites is growing as technologies for the CubeSat and associated hardware advance. An autonomous proximity operation was implemented to improve the performance of small satellites. The objectives of nanosatellite missions are divided into science-driven research and technological demonstrations. CanX-4 and -5 (Canadian Advanced Nanosatellite eXperiment-4 and -5) missions demonstrate precise formation-flying technology, such as drift recovery and station-keeping, formation-keeping, and reconfiguration with a dual-satellite [6,7]. GOMX-4 aims to demonstrate the formation-flying and inter-satellite link (ISL) technologies to prove the capabilities for future constellations [8,9]. The CubeSat proximity operation demonstration (CPOD) validates proximity operation technologies that consist of rendezvous, proximity operations, and docking (RPOD) with miniaturized technologies [10,11]. The objective of the Small-scale Magnetosphere and Ionosphere Plasma Experiment (SNIPE) mission is to explore the small-scale structure of terrestrial physical phenomena by collecting scientific temporal and spatial data [12,13,14].

1.2. Virtual Telescope

A virtual telescope (VT) composed of separate satellites, as opposed to a standard telescope, can not only change focal lengths but also demonstrate a longer focal length without manipulating the telescope’s structure. Depending on the type of VT, lenses and detectors are used in one pair [15,16] and the other is the (external) occulter and detector [17,18]. The detector of the latter is sometimes referred to as an optical instrument or a telescope. The separated satellites are inertially aligned in a “string-of-pearls” configuration relative to an inertial frame to create a VT; this is known as inertial alignment hold (IAH). In general astrophysical research, it is possible to demonstrate a relative distance of approximately 1000–72,000 km [15,17]. In solar science, X-ray images can be captured at 50–500 m [16,18].

1.2.1. Long-Baseline Telescope

To achieve longer focal lengths in the X-ray and gamma-ray regions of the electromagnetic spectrum, formation flying can be advantageous [19]. There have been numerous proposals for X-ray imaging missions, including XEUS (X-ray Evolving Universe Spectroscopy), Simbol-X, and FLIP-3 [16,20,21]. The XEUS was intended to place dual spacecraft in an L2 orbit, but it was absorbed into the Constellation-X project and transformed into a monolithic mission. The Simbol-X mission was supposed to operate in a highly elliptical orbit with a 4-day time span, but this plan was abandoned owing to budgetary constraints. For FLIP-3, the X-ray imaging lenses have a focal length of 100 m, which corresponds to the distance between the lens and the detector [16]. At the Sun–Earth L1 Lagrangian point, quasi-continuous aligned observation was attainable with a control accuracy of within 2 cm in the transverse direction and 25 cm in the axial direction. Milli-Arc-Second Structure IMager (MASSIM) is a mission concept that employs transmissive, refractive, and diffractive optics and aims for X-ray astrophysical measurements in space [15,19,22]. To achieve a milliarcsecond angular resolution, the MASSIM requires 1000 km of inter-spacecraft separation, within 25 mm of transverse displacement, and 10 km of axial offset. The formation will be located at L2 and will have 1–10 days to observe the celestial target. Fresnel and MAXIM were created as extreme versions of MASSIM, with the goal of achieving microarcsecond angular resolution [15,19,20,21,22,23,24]. With a 5–10 m-diameter lens and a focal length of ${10}^5$–${10}^6$ km, gamma-ray detection will examine the event horizon surrounding supermassive black holes.

1.2.2. Coronagraph with an External Occulter

The New World Observer (NWO) mission is a novel concept mission for exoplanet investigation involving two spacecraft [17,25,26]. External occulters are known to have superior scientific capabilities in internal coronagraphs and indirect methods. Starshade, which is equipped with a 50 m-diameter occulter, follows the telescope in the vicinity of the Earth–Moon–Sun L2 libration point [27,28]. The Starshade was held within 1 m of the line of sight of the telescope to the target star at a distance of 72 Mm. The alignment error is equivalent to 3 milliarcseconds, and the requirement for range control at 100 km is relatively low. Because the Telescope and Starshade orbit a perturbed heliocentric orbit at the unstable L2 point, the Starshade slews 1.4° per day owing to the dominant perturbations of the solar pressure and gravitational attraction. Consequently, the relative gravitational dynamics and drift are analyzed using the restricted N-body problem in the fields of the Sun and Earth’s gravitational Moon. The PROBA-3 (PRoject for On-Board Autonomy-3) mission, scheduled for launch in 2024, demonstrates and validates the solar coronagraph, a “large, virtually rigid structure” [18,29]. The two PROBA-3 spacecraft in formation produce the ASPIICS (Association of Spacecraft for Polarimetric and Imaging Investigation of the Corona of the Sun) coronagraph by utilizing the coronagraph spacecraft (CSC) with an optical instrument and the occulter spacecraft (OSC) with an occulting disk [30]. The orbit is designed to be a highly elliptical geosynchronous orbit with an inclination of approximately 59°, 600 km perigee altitude, and 60,530 km apogee altitude. The high eccentricity guarantees a six-hour-long observation at the apogee. The primary objective is to restrict the relative displacement error (RDE) to a few millimeters and the absolute attitude error (AAE) to a few arcseconds [31]. The formation is maintained by an H-infinity controller based on the state of unperturbed eccentric orbits, and it commands precision cold gas thrusters to compensate for six degrees of freedom of perturbations. This mission will lay the groundwork for autonomous formation-flying technologies to be used in future VTs [17,29]. By operating in a low-Earth orbit, a solar coronagraph can study solar eruptive events (SEEs), such as the release of coronal energy, sites where particles accelerate, and fundamental physical processes [32].

1.3. CANYVAL-C Mission

The CANYVAL project (CubeSat Astronomy by NASA and Yonsei using Virtual ALignment) was developed by Yonsei University and NASA and has the ultimate goal of demonstrating a VT by operating 1U and 2U CubeSats in low-Earth orbit. CANYVAL-X, the first phase of the project, was designed to implement the IAH technology experiment by constructing a VT configuration with a variable focal length [33,34]. As an advanced mission, CANYVAL-C aligns dual satellites at a fixed nominal distance such that the portrait of the corona is valid for executing a scientific mission [35].

1.4. Contribution

This paper describes the orbit control of the formation-flying technology developed for the CANYVAL-C mission. The mission operates the satellites in low-Earth orbit, maintaining a string-of-pearls configuration rather than a typical circular formation. Hence, the reference trajectory is defined by describing the observation target in a relative coordinate system, and the consecutive observation periods and alignment possibilities are evaluated based on this information. Because of the limited performance of cube satellites, control algorithms and complementary strategies have been developed and validated. The detector satellite is constructed as 1U, the minimum standard for cube spacecraft, whereas the occulter satellite has a propulsion system and an occulter in the 2U standard. A strategy was devised to minimize the degradation of mission performance due to differences in perturbations, and the mission success rate was verified to be 3 sigma or higher. Because the orbit control system was meant to be flight proven, Monte Carlo simulations were employed to study the effects of the other components of the GNC system, such as the navigation, attitude determination, attitude control, and propulsion system.

1.5. Brief Overview

Section 2 chronologically addresses the mission operations concept and the relative orbit control scenario, as well as the key requirements of this mission. The GNC architecture and a summary of the algorithms and requirements are concisely presented in Section 3. The relative orbit control algorithms and strategic approach for this mission are described and verified through software simulation in Section 4. Section 4.1, Section 4.2 and Section 4.3 cover the rendezvous, differential air drag control, and IAH, respectively. Section 5 presents the Monte Carlo simulation results for the entire scenario, and conclusions are summarized in Section 6.

2. Mission Description

2.1. Concept of the Mission

For the CANYVAL-C mission, the two CubeSats were in operation in a sun-synchronous orbit at an altitude of 525–550 km, an inclination of 97.6°, and a local time of ascending node (LTAN) of 11:00. The relative motion is described in the local-vertical and local-horizontal coordinate frame, which is a chief-centered rotating frame (Figure 1a). Because the relative displacement is less than tens of kilometers, it is represented by the Hill–Clohessy–Wiltshire equations (HCW), which are linearized equations of relative motion. The 1U CubeSat, the chief satellite, performs an attitude pointing toward the observation target, carrying an optical camera. The four nozzles of the 2U propulsion system can fire with two degrees-of-freedom perpendicular to the occulter (Figure 1b). The deployable occulter that is contained by the deputy satellite is 0.5U in size and has a cross-sectional area of $0.75\times 0.75$ m² when deployed.

2.2. Relative Orbit Control Scenario

Multiple formation flight techniques are required to achieve the ultimate goal of the CANYVAL-C mission and the IAH. The relative orbit control for this mission is chronologically shown in Figure 2. After the launch of the two satellites, during the Launch and Early Orbit Phase the satellites will be separated by thousands of kilometers owing to orbital perturbations during antenna deployment, hardware commissioning, angular velocity stabilization, data transmission and reception, and so on. When the relative distance between the two satellites is reduced to a few kilometers, their orbits can be described as a relative motion. Drift recovery and station-keeping (DRSK) regulates the satellite’s orbit using Gauss’ variational equation, enabling the deputy satellite to stop within 2 km of the relative distance from the chief satellite [36]. This study only deals with the relative orbit control algorithms immediately after DRSK.

The deputy satellite is then placed in a stable periodic relative orbit using rendezvous (RDV) within 0.5 T, where T is the orbital period, and the periodicity is maintained using differential air drag control (DADC) maneuvers. An IAH is performed if the satellite’s state is sufficient for the mission and the starting conditions for the mission mode are met. With the occulter stowed, the mission mode is executed once (commissioning-IAH) and the spacecraft is charged for five orbital periods while pointing toward the sun. In sun-pointing mode, the trajectory is uncontrolled, resulting in a drift from the perturbations. Consequently, the orbit control mode of (1) RDV and (2) DADC from Figure 2 is repeated ((5) RDV, (6) DADC). The deployment of the occulter increases the differences in perturbations; therefore, the occulter is deployed one day prior to the IAH to mitigate the drift.

2.3. Mission Requirement

The occulter must maintain a configuration with a nominal separation of 40 m from the observation target, the Sun, because the focal length of the optical camera is designed to be 40 m. Alignment must be performed during daylight, which has less than one orbital period. Therefore, the relative position, r, must be controlled to the desired trajectory, ${\mathit{r}}_d$, during the shaded daylight period (Figure 3). The control accuracy requirement is located within 3 m for each axis at the nominal position. This requirement was designed considering the size of the solar corona and the performance of the detector. Since the thrust system has only two degrees of freedom, it is impossible to fire the thrust in the direction of the occulter, making it impossible to control the relative distance. The requirements are divided into three groups: fine requirement, coarse requirement (1), and coarse requirement (2), as shown in

Table 1. IAH requirements for CANYVAL-C mission [35].

	Requirement
Nominal distance	40.0 m
Fine requirement	< 3.0 m for each axis
Coarse requirement (1)	${∥\mathit{r}-{\mathit{r}}_d∥}_2$ < 5.0 m	${\theta }_{align}$ < 7.5 degrees
Coarse requirement (2)	${∥\mathit{r}-{\mathit{r}}_d∥}_2$ < 10.0 m

. The most stringent requirement for fine alignment must be met to obtain an ideal depiction of the solar corona. Even if the relative distance offset is substantial, the solar corona can still be detected if the alignment angle is less than 7.5°. Coarse requirements (1) and (2) are defined as relative distance errors of 5 and 10 m, respectively, with an alignment angle under 7.5°. If these requirements are accomplished, the mission will be successful and the solar corona will be able to be obtained. The more detailed derivation of the mission requirement associated with the celestial object and payload is suggested in [35].

3. Solar Coronagraph Trajectory Design

In low-Earth orbit, the Sun is an inertial vector during a couple of orbital periods, so it shows a constant cross-track component. Therefore, the out-of-plane and in-plane motions of the Sun must be distinguished, and the less out-of-plane motion, the more efficient the staging mode. Assuming a sun-synchronous orbit, the changeable orbital elements are the semi-major axis and RAAN ($\Omega$, Right Ascension of the Ascending Node). The altitude was set as 550 km, so the influence of the RAAN should be determined. The RAAN, with respect to the Sun, can be expressed as the LTAN. Therefore, regardless of the position of the Sun, the closer the LTAN is to 00:00 or 12:00, the better the alignment.

The position of the Sun in the Earth-centered inertial frame (ECI frame) is expressed as an ecliptic plane. The vernal equinox, the intersections of the ecliptic plane, and the Earth’s equator plane define the x-axis of the ECI frame (Figure 4a). Because the maximum declination of the Sun for one year is 23.6°, the z-component of the Sun is always smaller than the magnitude of the vector projected onto the XY plane. The launch date of the CANYVAL-C mission was planned to be March, so the position of the Sun was set at the vernal equinox. Figure 4b shows the geometrical relationship of the sun-synchronous orbit according to the LTAN and Sun vector. If the LTAN differs by 12 h, the ascending and descending nodes appear on opposite sides of the same plane. According to the definition of the LTAN, the position of the ascending node with respect to the Sun ($\delta \Omega$) differs by 15° for every hour of the LTAN. Because the inclination is not exactly 90° for the actual sun-synchronous orbit, the angle between the orbital plane and the Sun ($\beta$) is calculated and presented in

Table 2. Beta angle and the cross-track ratio of the sunlight.

	LTAN 01	LTAN 03	LTAN 05	LTAN 07	LTAN 09	LTAN 11
$\delta \Omega$ (deg)	−15.00	−45.00	−75.00	−105.00	−135.00	−165.00
$\beta$ (deg)	−13.00	−42.65	−71.53	−105.14	−133.66	−163.28
$\sin \left(\beta \right)$	−0.23	−0.68	−0.95	−0.97	−0.72	−0.29
	LTAN 13	LTAN 15	LTAN 17	LTAN 19	LTAN 21	LTAN 23
$\delta \Omega$ (deg)	+15.00	+45.00	+75.00	+105.00	+135.00	+165.00
$\beta$ (deg)	+13.01	+42.66	+71.55	+105.11	+133.65	+163.27
$\sin \left(\beta \right)$	+0.23	+0.68	+0.95	+0.97	+0.72	+0.29

. An orbit with a large along-track component was the case for LTAN at 01, 11, 13, and 23. This indicates a difference of less than 3.5° between $\delta \Omega$ and $\beta$. $\sin \beta$ is the cross-track ratio of the Sun vector, ranging from 0.23 to 0.97.

The reference trajectory considering the position of the Sun in the LVLH frame is presented in Figure 5. From the left panel of Figure 5, the closer the LTAN is to 00 or 12, the smaller the cross-track component, and the closer it is to 06 and 18, the larger the cross-track component. According to the right panel in Figure 5, the cross-track component is constant, and the in-plane components change continuously.

Table 3. Daylight duration and daylight duration per orbital period at the vernal equinox.

	LTAN 01	LTAN 03	LTAN 05	LTAN 07	LTAN 09	LTAN 11
$T_{day}$ (min)	60.55	65.05	95.64	95.64	66.38	60.82
$T_{day}/T$ (%)	63.31	68.01	100.00	100.00	69.41	63.59
	LTAN 13	LTAN 15	LTAN 17	LTAN 19	LTAN 21	LTAN 23
$T_{day}$ (min)	60.58	65.12	95.64	95.64	66.32	60.78
$T_{day}/T$ (%)	63.34	68.08	100.00	100.00	69.34	63.55

shows the daylight duration ($T_{day}$) and ratio of the daylight duration to the orbital period ($T_{day}/T$). The orbital period was 95.64 min for every orbit, and the daylight period was approximately 63% for the LTAN 01, 11, 13, and 23 orbits and 100% for the LTAN 05, 07, 17, and 19 orbits. If the cross-track component of the LTAN is large, it shows a dawn–dusk orbit with long exposure to the Sun. As a result, there is no eclipse, which suggests many opportunities for alignment with the Sun. On the other hand, when the cross-track component is small, alignment in the along-track direction is easy, but the daylight duration is shortened, resulting in fewer opportunities for alignment.

4. GNC Subsystem Architecture

4.1. Navigation

The global navigation satellite system (GNSS)-based orbit estimation is performed autonomously onboard [37]. The orbital states processed from the GNSS receiver were used as the absolute orbit determination results and for the coarse relative navigation mode by simple subtraction. When a precise navigation result is required, it is determined using the differential global positioning system (DGPS) algorithm from the GNSS raw data. The requirements for the navigation subsystem are listed in

Table 4. Requirement of navigation subsystem [37].

Absolute Orbit Determination	Relative Orbit Determination
	Coarse Solution	Fine Solution
<10.0 m, 0.1 m/s	<20.0 m, 0.2 m/s	<1.0 m, 0.01 m/s
($3\sigma$, for each axis)	<($3\sigma$, for each axis)	($3\sigma$, for each axis)

.

4.2. Attitude Determination and Control

The attitude determination and control (ADCS) subsystem includes a fine sun sensor, photodiodes (coarse sun sensors), magnetometers, and gyro sensors, which estimate the attitude states as the quaternion and angular velocity using the extended Kalman filter (EKF) algorithm [38]. The deputy satellite, which requires precise attitude and orbit control, has two actuators: the reaction wheel assembly and the magnetorquer. The primary actuator is a reaction wheel that controls attitude pointing for charging, communication, and thrust firing. Because the reaction wheels generate torque through momentum conservation, the magnetorquer reduces the exceeded angular speed of the wheels. The requirements of the ADC subsystem are presented in

Table 5. Requirements of attitude determination and control subsystem [38].

Attitude Determination	Attitude Control
<1.5 deg, 0.15 deg/s	<3.0 deg, 0.5 deg/s
($3\sigma$, for each axis)	($3\sigma$, for each axis)

.

4.3. Guidance and Control

4.3.1. Propulsion System

The propulsion system used for the CANYVAL-C mission and equipped in the 2U was the Nanoprop of GomSpace. It has four nozzles at the center of each $xz$ and $yz$ plane and can fire in ${\widehat{b}}_x$ and ${\widehat{b}}_y$ directions of the body frame. The maximum magnitude of the thrust was 1.0 mN for each nozzle, and one or two nozzles were simultaneously fired towards the desired direction. Because the nozzles lie on the $xy$ plane of the body frame, the z-component of the control error is uncontrollable. Hence, for an orbit maneuver with three degrees of freedom, the attitude of the spacecraft must be adjusted simultaneously. The target thrust vector is projected onto the thrust plane, and the error angle between the target vector and the projected vector is calculated and controlled by the reaction wheel torque. According to the specifications of the reaction wheels, the maximum torque is 1.414 Nm and the spacecraft is controlled within 1.0°/s (

Table 6. Specification and configuration of the propulsion system.

Maximum Thrust	Thrust Resolution	Nozzle Configuration	Firing Accuracy ($3\mathit{\sigma }$)
1.0 mN	0.01 mN	$\pm {\widehat{b}}_x$, $\pm {\widehat{b}}_y$	<3.0 deg, 15%

).

4.3.2. Separation of Control and Estimation Sampling Intervals

The CubeSat was designed based on the 1U standard of 10 cm × 10 cm × 10 cm and has a very narrow gap between the sensor and other components. In such a configuration, the concurrent operation of any hardware increases sensor noise and degrades the estimation results. By separating the estimation and control intervals, the actuator should not be activated during attitude estimation. According to [38], when the attitude is estimated using the 2U satellite’s fine sun sensor the error converges within 2–3 s. The sampling rate of the orbit control input was set to 0.2 Hz and the attitude was estimated 2 s prior to sampling the thrust vector and then controlled for 3 s immediately after the sampling. The attitude observation data for the estimation and control torque were sampled at 1 Hz. When tracking the thrust vector for orbital maneuvering, the reference angular velocity is zero, and the current attitude does not change significantly if torque is not applied during the estimation interval. However, the relative velocity is not zero during relative orbit control; therefore, the zero-control input causes the deputy to drift away from the chief. Therefore, the thrust must be fired continuously, while the attitude is adjusted from its current state to the desired state.

5. Guidance and Control Algorithm

5.1. Rendezvous

The ideal goal of DRSK is to achieve zero relative velocity at a specific distance. However, orbit control errors result from the orbital prediction, ADC, and thrust errors when a real satellite is in operation. After DRSK is complete, it has a random relative distance and relative velocity within the specified range. To prevent the deputy from drifting away and settling in the projected circular orbit (PCO), it is necessary to perform rendezvous after station-keeping is complete.

If DRSK satisfies the requirements, the deputy is located within the relative position of 1 km ($3\sigma$) in the positive along-track direction and has a relative velocity of less than 0.12 m/s in a random direction. Because the magnitude of the relative velocity is obtained by multiplying the relative distance by the orbital angular velocity ($\omega =\sqrt{\mu /a^3}$), where $\mu$ is the gravitational standard constant and a is the semi-major axis, the reference trajectory is set as a PCO with a relative orbital radius of $1/\omega$ times the initial relative speed. Rendezvous is usually performed by fuel-optimal control based on the HCW equation. In this study, sliding mode control was additionally applied to control the errors from the nonlinear dynamics and disturbances. When operating a satellite, the relative orbit state is obtained from the DGPS algorithm, utilizing the GPS raw data from the chief. Therefore, a strategy is required to complement these errors in the guidance trajectory.

5.1.1. Fuel-Optimal Control Algorithm

According to [39,40], the analytic solution of the fuel optimization problem is derived using the calculus of variations and is described in the HCW equation. This analytical solution reduces the computational burden and explicitly obtains optimal solutions without solving the problem onboard. Consequently, the controlled reference orbital states and optimal thrust vector are derived as functions of the state vector and fundamental matrix. The general state-space equation is expressed as

$$\dot{\mathbf{X}}\left(t\right)=\mathbf{A}\left(t\right)\mathbf{X}\left(t\right)+\mathbf{B}\left(t\right)u\left(t\right),$$

(1)

where $\mathbf{X}\left(t\right)={\left[\begin{array}{cc}\mathit{r}\left(t\right)& \dot{\mathit{r}}\left(t\right)\end{array}\right]}^T$ is the state vector, $u\left(t\right)$ is the control input, $\mathbf{A}\left(t\right)$ is the state matrix, and $\mathbf{B}\left(t\right)$ is the input matrix. The state matrix $\mathbf{A}=\left[\begin{array}{cc}{\mathbf{0}}_{3\times 3}& {\mathbf{I}}_{3\times 3}\\ {\mathbf{A}}_1& {\mathbf{A}}_2\end{array}\right]$, ${\mathbf{A}}_1=\left[\begin{array}{ccc}3n^2& 0& 0\\ 0& 0& 0\\ 0& 0& -n^2\end{array}\right]$, and ${\mathbf{A}}_2=\left[\begin{array}{ccc}0& 2n& 0\\ -2n& 0& 0\\ 0& 0& 0\end{array}\right]$. The fundamental matrix associated with the HCW equation is well known from [41]:

$$\mathbf{\Phi }=\left[\begin{array}{c}{\mathbf{\Phi }}_{\mathbf{A}}\\ {\dot{\mathbf{\Phi }}}_{\mathbf{A}}\end{array}\right]=\left[\begin{array}{cccccc}4-3\cos nt& 0& 0& \frac{\sin nt}{n}& \frac{2(1-\cos nt)}{n}& 0\\ 6(\sin nt-nt)& 1& 0& \frac{2(\cos nt-1)}{n}& \frac{4\sin nt}{n}-3t& 0\\ 0& 0& \cos nt& 0& 0& \frac{\sin nt}{n}\\ 3n\sin nt& 0& 0& \cos nt& 2\sin nt& 0\\ 6n(\cos nt-1)& 0& 0& -2\sin nt& 4\cos nt-3& 0\\ 0& 0& -n\sin nt& 0& 0& \cos nt\end{array}\right],$$

(2)

where n denotes the mean motion.

For the fuel-optimal reconfiguration problem in circular orbits, the independent variable is time, where $b\left(t\right)=1$, $k=1$, and $\mathit{f}\left(t\right)=1$. Therefore, the cost function and analytical solution are derived as follows:

$$J=\frac{1}{2}{\int }_{t_0}^{t_{\mathit{f}}}{u}_{opt}^T{u}_{opt}d\tau .$$

(3)

Consequently, the explicit analytic solution is described as

$$\begin{array}{cc}\hfill \mathbf{S}& ={\int }_{t_0}^{t_{\mathit{f}}}{\mathbf{\Phi }}_{\mathit{A}}^T{\mathbf{\Phi }}_{\mathit{A}}d\tau \hfill \\ \hfill K& ={\mathbf{\Phi }}_{\mathit{f}}^{-1}{\mathbf{\xi }}_{\mathit{f}}-{\mathbf{\Phi }}_0^{-1}{\mathbf{\xi }}_0\hfill \\ \hfill u_{opt}& ={\mathbf{\Phi }}_{\mathit{A}}{\mathbf{S}}_{\mathit{f}}^{-1}\mathbf{C}K\hfill \\ \hfill J& =\frac{1}{2}{K}^T{\mathbf{C}}^T{\mathbf{S}}_{\mathit{f}}^{-1}\mathbf{C}K\hfill \\ \hfill =\left[\begin{array}{c}\mathit{r}\\ \dot{\mathit{r}}\end{array}\right]& =\mathbf{\Phi }({\mathbf{\Phi }}_0^{-1}{\mathbf{\xi }}_0+{\mathbf{C}}^{-1}\mathbf{S}{\mathbf{S}}_{\mathit{f}}^{-1}).\hfill \end{array}$$

(4)

5.1.2. Sliding Mode Control

The linear sliding variable is utilized to compensate for errors from nonlinear terms and perturbations and is designed as

$$\mathit{s}={\dot{\mathit{r}}}_e+\lambda {\mathit{r}}_e,$$

(5)

where the error vector is ${\mathit{r}}_e=\mathit{r}-{\mathit{r}}_d$ and its time derivative is ${\dot{\mathit{r}}}_e=\dot{\mathit{r}}-{\dot{\mathit{r}}}_d$.

If $\mathit{s}=0$, then the Lyapunov function $V=\frac{1}{2}{\mathit{r}}_e^T{\mathit{r}}_e$ derives the time derivative $\dot{V}=-2\lambda V$. This guarantees that the error converges to zero when the sliding variable is equal to zero.

The time derivative of the sliding variable is derived as

$$\begin{array}{cc}\hfill \dot{\mathit{s}}& ={\ddot{\mathit{r}}}_e+\lambda {\dot{\mathit{r}}}_e\hfill \\ \hfill & =\ddot{\mathit{r}}-{\ddot{\mathit{r}}}_d+\lambda {\dot{\mathit{r}}}_e\hfill \\ \hfill & =\mathit{f}+\frac{1}{m}(u+d)-{\ddot{\mathit{r}}}_d+\lambda {\dot{\mathit{r}}}_e\hfill \\ \hfill & =\mathit{\zeta }+\frac{1}{m}u,\hfill \end{array}$$

(6)

where $\mathit{\zeta }=\mathit{f}+\frac{1}{m}d-{\ddot{\mathit{r}}}_d+\lambda {\dot{\mathit{r}}}_e$ includes the nonlinear dynamics $\mathit{f}\left(t\right)$ and external disturbance term $d\left(t\right)$. $\mathit{\zeta }$ is assumed to have a known bound, D.

Let us define the sliding mode control law $u$ as

$$u_{SMC}\left(t\right)=-K\mathit{s}\left(t\right),$$

(7)

where K denotes constant control gain. Lyapunov theory was utilized to prove the stability of the dynamic system using the control input.

$$V=\frac{1}{2}{\mathit{s}}^T\mathit{s}$$

(8)

By differentiating the Lyapunov function (V), the time derivative of the Lyapunov function is derived as follows:

$$\dot{V}={\mathit{s}}^T\dot{\mathit{s}}.$$

(9)

Substituting Equations (6) and (7) into Equation (9) yields

$$\begin{array}{cc}\hfill \dot{V}& ={\ddot{\mathit{s}}}^T\dot{\mathit{s}}\hfill \\ \hfill & ={\ddot{\mathit{s}}}^T\left(\mathit{\zeta }-\frac{K}{m}\mathit{s}\right)\hfill \\ \hfill & <\ddot{∥\mathit{s}∥}\left(D-\frac{K}{m}∥\mathit{s}∥\right)\hfill \end{array}$$

(10)

where $∥·∥$ is a vector two-norm and D is the known bound of $\mathit{\zeta }$. Assuming that the sliding variable is in the region $∥\mathit{s}∥>\epsilon$, then

$$\begin{array}{c}\hfill \dot{V}<-\beta V^{1/2}\end{array}$$

(11)

where $\beta =\sqrt{2}·(\frac{K}{m}·\epsilon -D)$ is positive. Therefore, if the control gain is chosen as $K>\frac{Dm}{\epsilon }$, the sliding variable is bounded as $∥\mathit{s}∥\le \epsilon$.

5.1.3. Guidance Trajectory Update Strategy

The primary objective of rendezvous is to minimize the drift rate in the final state by generating a passively safe orbit. Therefore, a PCO with an initial position and velocity was generated, and the final state of the PCO was set as the target state. During half of the orbital period for rendezvous, the initial time and state are updated and a new optimal trajectory is generated for every cycle (120 s), which is the ISL period. As a result, the periodic guidance update strategy not only compensates for the error due to the nonlinearity of dynamics but also prevents the initial estimated error that is contained in the optimal trajectory from gradually increasing as it propagates. The formulation for PCO was obtained from [42]:

$$\left\{\begin{array}{c}x\left(t\right)=\frac{d_1}{2}\sin (nt+\alpha )\hfill \\ y\left(t\right)=d_1\cos (nt+\alpha )+d_3\hfill \\ z\left(t\right)=d_2\sin (nt+\beta )\hfill \end{array}\right..$$

(12)

The guidance trajectory of the sampling frequency (${\mathit{f}}_{S1}$) was chosen as 2 to 10 min (Figure 6). If it is set to 2 min, the final drift rate may decrease, but it requires a larger thrust owing to its shortened maneuver time. In this study, ${\mathit{f}}_{S1}$ is set as 0.0083 Hz (once per 120 s) because the maximum thrust is smaller than the saturation limit, and the fuel budget is affordable. During orbit control, the target thrust vector should not change frequently; otherwise, the attitude control performance should be degraded. Hence, the thrust vector sampling frequency (${\mathit{f}}_{S2}$) is set as 0.25 Hz (once per four seconds) and fired consistently.

Figure 7 shows the rendezvous simulation results for achieving the PCO. The dynamic models were set as JGM3 70 × 70 for gravity, the exponential model for atmospheric drag, and the spherical model for solar radiation pressure. The differential equation was integrated using the ode45 integrator in MATLAB with a step size of 1 s. The simulation settings were identical to those used in this study. In the final controlled state, the error is approximately 0.1 m in the radial and cross-track directions and 1.0 m in the along-track direction. The optimal control input shows discontinuity every 120 s, and the magnitude increases after half the maneuver time. This is because the remaining feasible time is reduced and can be eliminated when ${\mathit{f}}_{S1}$ is decreased. While the optimal control input changes over time, the sliding mode control input has a similar form to noise under 0.02 mN. Initially, it barely required a cross-track maneuver; thus, the optimal control showed a lower level than the sliding mode control.

5.2. Differential Air Drag Control

After rendezvous, satellites in periodic orbit drift by perturbation. Therefore, it is necessary to control the formation of satellites to maintain stability. The difference in atmospheric drag between the two satellites allows the formation to be maintained for months without thrust consumption. In this study, differential atmospheric drag control was performed by controlling the radial and along-track positions that had the greatest impact on drift.

Averaged In-Plane and Drift Terms

Factors that affect differential air drag include the structural properties of the spacecraft, atmospheric density, solar flux, geomagnetic activity, and operational altitude [43]. The air drag is described as

$${\mathit{a}}_{drag}=-\frac{1}{2}\frac{C_D\mathit{A}}{m}\rho V_{rel}^2{\widehat{\mathit{V}}}_{rel},$$

(13)

where $C_D$ is the drag coefficient, A is the effective cross-sectional area, m is the mass of the spacecraft, $\rho$ is the atmospheric density, and $\mathit{V}=V_{rel}{\widehat{\mathit{V}}}_{rel}$ is the velocity relative to the atmosphere. The differential air drag that occurs when these factors are different is described as

$$\ddot{y}\left(t\right)=\frac{3}{2}{V}_{rel}^2\rho C_D(\frac{{\mathit{A}}_{deputy}}{m_{deputy}}-\frac{{\mathit{A}}_{chief}}{m_{chief}}).$$

(14)

Atmospheric density is affected by solar flux, geomagnetic activity, and altitude. The intensity of the solar flux is expressed as $F10.7$, and the geomagnetic activity is the $K_p$ or ${\mathit{A}}_p$ index. According to [43], the along-track separation increases nonlinearly as $F10.7$ and ${\mathit{A}}_p$ increase. These values change under the influence of solar activity, with a solar cycle of approximately 11 years, and show a similar level within a few months. Therefore, the air density was calculated using Equation (15) by obtaining observation data for a few months before the operation period from Celestrak.

$$\left\{\begin{array}{c}T=900+2.5(F10.7-70)+1.5{\mathit{A}}_p\hfill \\ m=27-0.012(h-200)\hfill \\ H=T/m\hfill \\ \rho =6\times {10}^{-10}exp(-(h-175)/H),\hfill \end{array}\right.$$

(15)

where T (Kelvin) is a fixed exospheric temperature, m is the actual variation in molecular mass, and H (km) is a variable scale height to describe the variation of the molecular mass, with the altitude h being 180 to 500 km.

Because the air drag depends on the velocity with respect to the atmosphere, the DADC strategy was established by analyzing the radial and along-track direction components of the HCW equation. The average value during one orbital cycle of the radial and along-track components was calculated using an averaging operator to remove the oscillatory motion as follows:

$$\begin{array}{cc}\hfill \overline{x}& =〈x\left(t\right)〉\triangleq \frac{1}{T}{\int }_0^{T}x\left(\tau \right)d\tau \hfill \\ \hfill & =\frac{1}{T}{\int }_0^T\left[4x\left(0\right)+\frac{2\dot{y}\left(0\right)}{n}\right]+\frac{\dot{x}\left(0\right)}{n}\sin n\tau -\left[3x\left(0\right)+\frac{2\dot{y}\left(0\right)}{n}\right]\cos n\tau d\tau \hfill \\ \hfill & =4x\left(0\right)+\frac{2\dot{y}\left(0\right)}{n},\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill \overline{y}& =〈y\left(t\right)〉\triangleq \frac{1}{T}{\int }_0^{T}y\left(\tau \right)d\tau \hfill \\ \hfill & =\frac{1}{T}{\int }_0^T-\left[6nx\left(0\right)+3\dot{y}\left(0\right)\right]\tau +\left[y\left(0\right)-\frac{2\dot{x}\left(0\right)}{n}\right]+\left[6x\left(0\right)+\frac{4\dot{y}\left(0\right)}{n}\right]\sin n\tau +\frac{2x\left(0\right)}{n}\cos n\tau d\tau \hfill \\ \hfill & =-\frac{1}{2}(6nx\left(0\right)+3\dot{y}\left(0\right))T+\left[y\left(0\right)-\frac{2\dot{x}\left(0\right)}{n}\right]\hfill \\ \hfill & ={\overline{y}}_1+{\overline{y}}_2.\hfill \end{array}$$

(17)

Periodic functions such as sine and cosine functions are removed, and as a result the remaining terms are only the radial and along-track initial distances. The transverse component has a drift term that increases linearly with time. The control drag area ($\delta \mathit{A}$) was determined by multiplying the three values by the appropriate gains.

$$\delta \mathit{A}\left(t\right)=K_1{\overline{y}}_1+K_2{\overline{y}}_2+K_3\overline{x},$$

(18)

where $K_{1,2,3}$ denotes the control gains.

For radial and along-track separation, the navigation error is generally greater than the nominal value, and the control state easily becomes noisy, which makes $K_2$ and $K_3$ small. In this control algorithm, it is only possible to switch the attitude orientation once every two to three days, owing to the limited power. For formation flying, the altitudes of the two spacecraft are slightly different, so the ballistic coefficient ($C_B$) has the most significant impact. The ballistic coefficient is determined by the drag coefficient and the area-to-mass ratio.

$$C_B=\frac{C_D\mathit{A}}{m}$$

(19)

The attitude of the deputy is controlled to the maximum area if the summation of the equality area and the control area is larger than the maximum, otherwise it is maintained as the minimum.

$${\mathit{A}}_{ct\mathit{r}l}={\mathit{A}}_{Eq}+\delta \mathit{A},$$

(20)

where ${\mathit{A}}_{Eq}=\frac{{\mathit{A}}_{chie\mathit{f}}}{m_{chie\mathit{f}}}·m_{deputy}$ is the area of the deputy with the same area-to-mass ratio as that of the chief.

When DADC was performed at a distance of approximately 1.4 km with an initial drift rate of 0.5 cm/s, the relative trajectory for the five orbital periods was expressed as shown in Figure 8. Because it is initially positioned in a positive along-track direction, it has a negative drift rate by controlling its attitude to have a minimum effective cross-sectional area (top panel). Owing to the differential air drag, the direction of the relative velocity of the satellite changes from negative along-track to positive when the control drag area reaches its maximum effective cross-sectional area. In the negative along-track direction, the direction of velocity changes at approximately 2.8 km, and differential gravity occurs due to relative distance, thereby distorting the initial shape of the orbit (second panel). In the direction of the positive along-track, the distance of the orbit is the greatest (third panel), as it is separated by up to 14 km, and the differential gravity is the greatest. Because DADC is a passive control, the maximum distance at zero relative velocity gradually increases over time as the control error caused by perturbations accumulates (second and bottom panels).

5.3. Inertial Alignment Hold

The IAH is an orbit maneuver for retaining a satellite’s relative position constant with respect to an inertial frame. The position of the 2U satellite carrying the occulter and thruster is constantly adjusted against the Sun so that the 1U spacecraft carrying the camera can take photos. Owing to the limited battery power and operational time, the 1U–2U formation is reconfigured in the staging mode and aligned with the Sun in the alignment mode. In the sun-synchronous orbit, the position of the Sun is described in the LVLH frame to designate a reference trajectory, and a control strategy for the staging mode and alignment mode is constructed from this.

5.3.1. Mode Separation in Inertial Alignment Hold

Because DADC is performed before IAH, the relative orbit cannot perfectly maintain its periodicity by perturbation. As the two satellites move toward and away, the navigation error makes the initial states of the mission mode uncertain. In addition, an eclipse should be avoided for alignment with respect to the Sun. This uncertainty inspires the consideration of an efficient strategy rather than immediate alignment. The strategy is to conform to a 1U–2U “string-of-pearls” configuration of 40 m (staging mode) and align with respect to the Sun right after it gets out of the eclipse (alignment mode). The relative motion with control acceleration is represented by [44]:

$$\left\{\begin{array}{c}\ddot{x}\left(t\right)=3n^{2}x+2n\dot{y}+u_x/m\hfill \\ \ddot{y}\left(t\right)=-2n\dot{x}+u_y/m\hfill \\ \ddot{z}\left(t\right)=-n^{2}x+u_z/m,\hfill \end{array}\right.$$

(21)

where x, y, z are the relative states in radial, along-track, and cross-track components in the LVLH frame, n is the mean motion, and $u={\left[\begin{array}{ccc}{u}_x& u_y& u_z\end{array}\right]}^T$ denotes the control input.

To maintain the relative configuration without control inputs during the staging mode, the radial and cross-track positions and the radial and along-track velocities should be zero ($x=z=0$ and $\dot{x}=\dot{y}=0$). In addition, if the cross-track velocity exists, its position will eventually increase, so z must be zero. Consequently, if only the position in the along-track direction exists it is possible to maintain a stable formation without consuming thrust. Therefore, it is recommended to maintain 40 m separation in the positive along-track direction. Because the deputy is aligned in the along-track direction during the staging mode, the closer the Sun is to the along-track direction in the alignment mode, the better the alignment will be performed. Consequently, the strategy of two-submode operation is more advantageous when the Sun’s cross-track component is small.

5.3.2. Linear Quadratic Regulator Algorithm

A linear quadratic regulator (LQR) is used to regulate the relative states in a constant position in the positive along-track direction. The cost function is defined as

$$J=\frac{1}{2}{\int }_0^{\infty }\left({\mathbf{X}}_e^T\mathbf{Q}{\mathbf{X}}_e+u_{LQR}^T\mathbf{R}{u}_{LQR}+2{\mathbf{X}}_e^T\mathbf{N}{u}_{LQR}\right)dt,$$

(22)

where ${\mathbf{X}}_e=\mathbf{X}-{\mathbf{X}}_d$ is the relative control error; $\mathbf{X}$ and ${\mathbf{X}}_d$ are the relative and desired states, respectively; and $u_{LQR}$ is the LQR control input.

The linearized system is described as

$${\dot{\mathbf{X}}}_e=\mathit{A}{\mathbf{X}}_e+\mathbf{B}{u}_{LQR},$$

(23)

where the state and input matrices are obtained from Equation (1).

The Riccati equation and the optimal gain, K, are obtained by

$${\mathit{A}}^T\mathbf{S}+\mathbf{S}\mathit{A}-(\mathbf{S}\mathbf{B}+\mathbf{N}){\mathbf{R}}^{-1}({\mathbf{B}}^T\mathbf{S}+{\mathbf{N}}^T)+\mathbf{Q}=\mathbf{0},$$

(24)

$$K={\mathbf{R}}^{-1}({\mathbf{B}}^T\mathbf{S}+{\mathbf{N}}^T).$$

(25)

Finally, the control input of LQR is derived as

$$u_{LQR}=-K{\mathbf{X}}_e=-{\mathbf{R}}^{-1}({\mathbf{B}}^T\mathbf{S}+{\mathbf{N}}^T){\mathbf{X}}_e,$$

(26)

where the matrix $\mathbf{Q}$, $\mathbf{R}$, and $\mathbf{N}$ are user-specified values.

The errors are described in the block diagram of the control algorithm in Figure 9. The relative navigation error is included in the state input and affects the guidance trajectory, control input, and target quaternion of the tracking problem. The pointing error of the ADCS system is generated by the current estimation attitude and degrades the performance at thrust firing and payload alignment. The thrust errors are modeled as the misalignment and magnitude error of the thrust vector and are assumed to be applied at the center of mass of the satellite.

6. Monte Carlo Simulation

6.1. Rendezvous

Figure 10a shows that the drift rate has a Gaussian distribution with an average of 0.26 and three-sigma of 4.61 cm/s after the rendezvous. The final drift rate ($\Delta {\dot{\lambda }}_{\mathit{f}}$) is decreased compared to the initial drift rate ($\Delta {\dot{\lambda }}_i$), with a probability of approximately 90% (Figure 10b, top). In the remaining cases, the thrusts were saturated and the final drift rate was not sufficiently reduced. The final drift rate and drift distance exhibit a linear relationship (Figure 10b, bottom). Finally, the maximum drift distance is 270 m, which is sufficient for performing the DADC.

6.2. Differential Air Drag Control

After rendezvous, the DADC is executed in the normal mode for proximity operations without propellant and power consumption. As described in Section 5.3.1, the IAH positions the deputy in the positive along-track direction (staging mode) and initiates alignment with the Sun as soon as the daylight interval begins (alignment mode).

Table 7. Starting conditions for C-IAH.

Mission Mode Starting Conditions
(1) Along-track distance ($\left|y\right|$)	60 < $\left|y\right|$ < 70 (m)
(2) Staging duration ($T_{stg}$)	40 < $T_{stg}$ < 60 (min)

presents the conditions for the effective commencement of the IAH. Because the entrance into mission mode is performed autonomously onboard, the coarse solution of relative navigation must fulfill condition (1) in Table 7. The staging duration, condition (2) of Table 7, is defined as the time between entering mission mode and sunrise, and the relative orbit control error is stabilized only when a minimum of 40 min is guaranteed.

The DADC was simulated with 100 samples over a period of two months to assess the number and duration of mission successes (Figure 11). The total number of mission chances must be at least one to complete the mission scenario. As a consequence of the simulation, all the cases had at least two mission opportunities, with the majority having six. When the mission chance has a maximum value of 31, the control performance is outstanding, even when the error is considered, resulting in an extensive period of satisfying the condition (1) given in Table 7. Owing to condition (2) in Table 7, the mission chance appears discontinuously and regularly over several orbital periods during this time. The duration of the mission chance is determined by the relative orbit phase angle, which affects the amount of time before sunrise, and the drift rate, which is the amount of time it takes for the along-track position to change by 10 m and shows an average performance of 3.69 min.

6.3. Commissioning Inertial Alignment Hold

As a consequence of the Monte Carlo simulation, the random initial states that satisfy Table 7 converge to the desired trajectory within the limited operation time in the staging mode (Figure 12a). The relative error states in the radial and cross-track directions are controlled in the alignment mode, and the states cannot be perfectly controlled because of the restricted degrees of freedom of the thruster (Figure 12b). One sample in Figure 12c shows the deviation in the alignment angle, and it is the failure case of the coarse requirement (2) in

Table 8. Percentage of the mission succeeded cases and the minimum and maximum aligned time for fine requirement, coarse requirement (1), and coarse requirement (2) during C-IAH.

		Fine Req.	Coarse Req. (1)	Coarse Req. (2)
Mission success		75%	92%	99%
Aligned time (min)	min.	0.08	0.12	2.15
	Max.	12.87	16.65	16.68

. Comparing the alignment times of the three requirements (Table 1), the alignment time and success probability are the greatest for the coarse requirement (2), coarse requirement (1), and fine requirement (Figure 13 and Table 8). This indicates a 99% chance of success based on coarse requirement (2). The probability of success is 99% based on coarse requirement (2), whereas the probability of success of the fine requirement is remarkably low at 75%. This is because this satellite has a two-degrees-of-freedom propulsion system; therefore, if it orients toward the Sun in the alignment mode the alignment angle is well controlled, but the relative distance is difficult to control. The average impulse consumption is 3.41 Ns, with an average of 2.38 Ns in staging mode and 1.02 Ns in alignment mode (

Table 9. Impulse consumption for the staging and alignment mode and total impulse during C-IAH.

		Staging Mode	Alignment Mode	Total
Impulse (Ns)	Avg. $\pm 3\sigma$	2.38 ± 2.68	1.02 ± 0.50	3.41 ± 2.74
	Max.	4.79	1.73	5.68

). Because the staging mode is 40 to 60 min and the alignment mode is 20 min, it exhibits a comparable level of impulse per time.

6.4. Inertial Alignment Hold

According to the relative orbit control scenario, after the first mission mode (commissioning-IAH) with the occulter stowed, the deputy must orient toward the Sun for five orbital periods to ensure a sufficient power budget. Subsequently, the rendezvous and DADC are executed sequentially with the occulter deployed.

Table 10. Structural properties of 1U and 2U [35].

Area-to-Mass Ratio (m2/kg)	Minimum	Maximum (Stowed)	Maximum (Deployed)
2U	0.9578 $\times {10}^{-2}$	3.2574 $\times {10}^{-2}$	23.7342 $\times {10}^{-2}$
1U	1.6719 $\times {10}^{-2}$
Differential area-to-mass ratio (m2/kg)	Minimum	Maximum (Stowed)	Maximum (Deployed)
2U-1U	−0.7141 $\times {10}^{-2}$	+1.5855 $\times {10}^{-2}$	+22.0623 $\times {10}^{-2}$

lists the structural properties based on deployment. After deployment, the maximum cross-sectional area increased by approximately sevenfold, whereas the differential air drag increased by approximately fourteenfold. To compensate for this, the deputy should deploy the occulter and enter mission mode in the along-track direction to maintain the minimum cross-sectional area (

Table 11. Starting conditions for the IAH.

Mission Mode Starting Conditions
(1) Along-track distance (y)	60 < y < 140 (m)
(2) Staging duration ($T_{stg}$)	40 < $T_{stg}$ < 60 (min)

(1)). The staging duration is the same as that in Table 7.

As a result of the Monte Carlo simulation, the final drift rate of rendezvous was 0.06 m/s or less, even after C-IAH and Sun-pointing mode. The number and duration of mission chances in the DADC for two months are shown in Figure 14. At least one mission chance is guaranteed, but because it enters the mission mode only in the positive direction it has a much smaller number of mission chances compared to Figure 11. However, the duration of the mission chance did not decrease significantly because the range of the along-track distance increased from 10 to 80 m.

Compared to C-IAH, IAH shows larger control errors in the staging mode, and some control errors even diverge in the alignment mode (Figure 15). Similar to Figure 12c, the six samples of the deviated alignment angle in Figure 15c are the failure cases of coarse requirement (2) in

Table 12. Percentage of the mission succeeded cases and the minimum and maximum aligned time for fine requirement, coarse requirement (1), and coarse requirement (2) during IAH.

		Fine Req.	Coarse Req. (1)	Coarse Req. (2)
Mission success		64%	87%	94%
Aligned time (min)	min.	0.03	0.07	0.85
	Max.	10.77	16.57	16.78

. Consequently, it has a lower probability of success for all the fine requirement and coarse requirements (1) and (2) (Figure 16, Table 12). As the effective cross-sectional area increases, the air drag applied to the satellite increases; thus, both the orbital perturbation and the attitude disturbance increase. The orbital perturbation by air drag only has a velocity direction and is controllable by thrust. However, the air drag disturbance torque continuously accumulates and is difficult to compensate for with the MTQ, which degrades the attitude control performance of the staging and alignment modes, resulting in a low probability of success. Because attitude control is necessary for thrust firing, thrust consumption also increases (

Table 13. Impulse consumption for the staging and alignment mode and total impulse during IAH.

		Staging Mode	Alignment Mode	Total
Impulse (Ns)	Avg. $\pm 3\sigma$	2.83 ± 3.59	1.07 ± 0.63	3.90 ± 3.87
	Max.	7.08	1.82	8.90

). However, because the success probabilities of C-IAH and IAH are 99 percent and 94 percent, respectively, the mission has a 99.94% chance of success and conforms to the 3-sigma rule.

6.5. Fuel Budget

The total impulse of the propulsion system was 40 Ns, and the propellant budget was calculated to be 20% of the margin (

Table 14. Fuel budget for the relative orbit control scenario in CANYVAL-C mission.

	Impulse (Ns)	Attempts (#)	Accumulated Impulse (Ns)	Accumulated Impulse with 20% Margin (Ns)
Drift Recovery *	0.1876	1	0.1876	0.1876
Station-keeping *	0.1796	1	0.1796	0.1796
Rendezvous	0.9211	2	1.8422	2.2106
Inertial Alignment Hold	3.8976	2	7.7952	9.3542
Total Impulse	-	-	10.0046	12.0055

* Ref. [35].

). This shows that the fuel is sufficient to perform the entire mission, even within the margin. The mass of the fuel is only a few grams, which is very small relative to the mass of the spacecraft, so the changes in mass due to thrust consumption are ignored.

7. Conclusions

In this study, relative orbit control algorithms and demonstration scenarios for a VT with CubeSats in low-Earth orbit are established and validated. A solar coronagraph is a VT that photographs the solar corona by aligning two satellites with respect to the inertial coordinate frame. When a satellite is operated in a sun-synchronous orbit, the position of the Sun relative to the satellite orbit is constant. Consequently, the reference trajectory was defined as a circular trajectory with a constant cross-track offset on a sphere of radius equal to the focal length. CANYVAL-C is a CubeSat mission that demonstrates the use of an IAH to implement a solar coronagraph. To accomplish this mission, relative orbit control algorithms such as rendezvous, differential air drag control, and inertial alignment hold have been developed, as well as supplementary strategies to compensate for the limited performance of cube satellites and the relative orbit control scenarios. To validate this research, the subsystem errors and operational constraints of the cube satellite were considered in Monte Carlo simulations. Consequently, it was verified that the mission had a 3-sigma probability of success.