PicSat's Enduring Legacy. Probing the Flight of a Small Astronomical Satellite

The time is ripe for the world science community to consider that small satellites may play important roles, specially in astrophysics, planetary science, heliophysics, and Earth science (Millan et al. 2019). Indeed, the last Decadal Surveys lead by the U.S. National Academies and published by NASA ⁴ recommended increasing the level of investments in small satellite programs, mainly in the above mentioned areas. ESA has also had many initiatives concerning small satellites in the last years. ⁵

CubeSats do not replace larger missions. They can be used advantageously to supply complementary measurements and explore targeted science goals. An important characteristic concerning the small satellites is the "fly-learn-refly" procedure, which is often adopted for when the budget allows it. A first flight model is launched and, if any serious problem appears, a second model is launched with the proper modifications. This approach increases significantly the number of successful missions. CubeSats thus bring versatility to space exploration. In astrophysics, in spite of the small dimensions which limit optical instrument aperture, long-term monitoring can be quite useful for stellar variability and exoplanetary studies. These were the main goals of PicSat.

Constellations of CubeSats with a variety of purposes, including long-term Earth monitoring and solar and space physics studies (e.g., Liddle et al. 2020, and references therein) have been launched during the last years (Kulu 2020).

Attitude and Orbit Determination and Control System is responsible for estimating and controlling nanosatellite attitude and orbit. It is divided into ADCS itself, and Orbit Determination and Control System. These subsystems are subdivided into more specific ones: Attitude Determination System (ADS), Attitude Control System (ACS), Orbit Determination System (ODS), and Orbit Control System (OCS). Such a division separates measurements and state estimation algorithm from control law.

Since each CubeSat unit weights no more than 1.33 kg, technology miniaturization is required, orbit control is a challenge, and three-axis stabilization took time to become common place, while high accuracy remains unusual. A survey of 426 CubeSats launched until 2016 January (Polat et al. 2016) found that 56% of these CubeSats used magnetic actuators, 44% used reaction wheels, 13% used passive magnetic control, and 4% who did not use any type of control. Only 44%, the ones using reaction wheels, could implement three-axis control. Most of the CubeSats did without star trackers, limiting pointing accuracy to the order of degrees.

Except for MarCO-A and MarCO-B, all CubeSats were launched in LEO, of which 76% were in the 350–700 km range (Polat et al. 2016; Villela et al. 2019). As pointed out in Section 1, the development of high-accuracy ADCS requires simulators. Dynamic simulation of orbit and attitude allows assessing the orientation, trajectory and environmental conditions to which the spacecraft will be subject over time. For Sebestyen et al. (2018), in LEO, high-fidelity mathematical models must contain at least: (i) rigid body dynamics; (ii) orbital mechanics model; (iii) high-fidelity magnetic field model; (iv) high-fidelity atmospheric density model; (v) torque due to atmospheric drag; (vi) torque due to the pressure of solar radiation; (vii) torque due to the gravity gradient; (viii) torque due to the residual magnetic dipole; and (ix) position of the Sun.

Several academic simulators, when compared to the expected characteristics of a high-fidelity mathematical simulator, are deficient in at least one of the above items. This occurs both in simulators for specific applications, such as ACS and/or ADS for a particular nanosatellite Krogh & Schreder 2002; Gießelmann 2006; Francois-Lavet 2010; Jensen & Vinther 2010; Holst 2014; van Vuuren 2015; Rondão 2016; Thomsen & Nielsen 2016; Avanzini et al. 2019; Jonsson 2019), and in simulators with more comprehensive proposals to serve as a development platform (Menges et al. 1998; Lovera 2003; Nasirian et al. 2006; Naqvi & Raza 2007; Triharjanto et al. 2015; Habibkhah et al. 2017; Alshamy et al. 2019; Liu et al.2019).

There exist open source simulators and toolboxes such as ODTBX, developed by NASA, which deals only with orbital dynamics; ⁶ KPS available in C ++, MATLAB and Python, which has disturbance torques limitations (Omar 2017); PROPAT available as a MATLAB toolbox, also with disturbances limitations (Carrara 2015); SNAP available in Simulink and MATLAB, which has only attitude dynamics and limited disturbance modeling (Rawashdeh 2019); and OPEN-SESSAME developed in C ++ (Turner 2003).

There exist proprietary software Systems Tool Kit and Orbit Determination Tool Kit simulators from AGI, for either standalone use or integrated with MATLAB; ⁷ and Spacecraft Control Toolbox from Princeton Satellite System, available as a MATLAB package. ⁸ Given that the mathematical models of these simulators are not available, it is not possible to compare them to high-fidelity simulators criteria. However, all of them report being high-fidelity mathematical simulators.

Most of the cited simulators have no in-orbit data validation (Menges et al. 1998; Lovera 2003; Turner 2003; Nasirian et al. 2006; Carrara 2015; Triharjanto et al. 2015; Habibkhah et al. 2017; Omar 2017; Rawashdeh 2019). This is mainly due to the difficulty of obtaining real data. Alternatively, some studies compare their own simulation results with other academic (Liu et al. 2019) or commercial (Naqvi & Raza 2007; Alshamy et al. 2019) simulators.

In Sternberg et al. (2018), a team from NASA Jet Propulsion Laboratory uses flight data from the CubeSat MicroMAS-1 to validate their simulator. The same can be seen with CubeSats UWE-3 (Bangert et al. 2014), and ESTCube-1 (Sato et al. 2017). Some works are dedicated to validating the performance of one of their subsystems using in-flight data, which is the case for nanosatellites RAX-1 and RAX-2 (Slavinskis et al. 2016) and ESTCube-1 (Springmann & Cutler 2014). Finally, some other works use in-orbit data for an overall performance analysis (Taraba et al. 2009; Scholz et al. 2010; Xiang et al. 2012; Dechao et al. 2014; Sarda et al. 2014; Mason et al. 2017; Pong 2018).

ADCS simulators validated with in-flight data are altogether rare in the literature, showing the importance of using PicSat in-orbit data for this purpose, as we present in this paper.

The retained in-house payload final design was based on an innovative technique of space-based interferometry: photometry was acquired by means of an optical system composed of two mirrors and a single-mode optical fiber, and a detector system: a solid-state photodetector Single-Photon Avalanche Diode (SPAD). Figure 2 shows the payload schematically.

Zoom In Zoom Out Reset image size

The planetary transit detection method required high-precision photometry, which in turn required high-accuracy pointing and SPAD thermal stability. The mission required pointing the satellite toward Beta Pictoris with about 1'' rms accuracy in the instrument boresight frame. Such accuracy could not be achieved using only COTS ADCS components. Coping with these requirements necessitated another innovation: the satellite fine pointing control was divided into two stages. The first one was a three-axis stabilization COTS ADCS using star tracker, reaction wheels and gyrometers at platform level; and the second one an in-house developed active correction of the fiber position using piezoelectric actuators at payload level. The piezoelectric stage was mounted around the fiber to control its position regarding the target star and controlled in closed loop via a microprocessed system, enabling arc second accuracy on star pointing.

An embedded thermal control allowed SPAD noise to operate under 40 ppm hr⁻¹, adequate for the overall photometric error budget (110 ppm hr⁻¹) (Nowak et al. 2017).

ADCS and piezo-controlled fiber formed a two-stage position control system. ADCS controlled satellite angular position up to 30'' accuracy, and piezos controlled the linear fiber position through the observed target starlight, allowing 1'' final accuracy.

Figure 3 presents a simplified control diagram of the system, indicating the relevant ACS and ADS components. Desired and actual states can be either angular position or angular velocity depending on the platform control mode (de-tumbling or target pointing mode). The ADCS (model iADCS 100) was developed by Hyperion Technologies. Table 1 lists its actuators and Table 2 its sensors, according to its datasheet.

Zoom In Zoom Out Reset image size

As iADCS-100 is a proprietary device, we do not have access to detailed information about its algorithms. Nevertheless, by the time of the mission (iADCS-100 software version 2.00), the target pointing mode was probably running a state estimator based on a linear Kalman filter and three Proportional-Integral-Derivative controllers (one for each axis) with anti-Windup (for the Integral term).

The ADCS residual pointing errors (30'' accuracy) are the input of the Piezo-Stage Control System, whose objective is to correct them up to 1'' accuracy. There is no feed feedback implemented between them. To do so, the fiber is mounted on a two-axis piezo actuator and positioned in the focal plane of the telescope. The piezo actuator is based on space-qualified components and is made by CEDRAT Technologies. Each axis is controlled independently. They are also equipped with strain gauges to measure the piezo excursion. The total range of each piezo is about 430 μみゅーm. The fiber position is modulated in the focal plane. After a first search to find the star in the field of view, multiple photometric measurements acquired at different positions along the modulation pattern are combined in an extended Kalman filter to derive a real-time estimate of the star position. The loop is closed using the best estimation of the star position as a reference position to the fiber (Nowak 2019).

Figure 4 presents platform and payload components and their locations inside the 3U envelope. The upper cube unit comprises the power subsystem (EPS), the communication subsystems (TRXVU and HDRM), the VHF antennas and the on-board computer (OBC). The intermediate unit comprises the UHF antennas, the ADCS and the electronic payload, including the photodiode. Optomechanical payloads (Piezo-stage and Telescope) and the star tracker are accommodated in the lower unit. Further details can be found in Nowak et al. (2018).

Zoom In Zoom Out Reset image size

Payload subsystems were in-house developments and platform components were all COTS.

PicSat was launched on 2018 January 12, by the PSLV-C40 (Polar Satellite Vehicle Launcher), operated by the Indian Space Research Organization. It was placed at a 505 km Sun-synchronous polar orbit (orbital period of about 95 minutes).

Two major failures prevented the mission from being fully successful: one in the ADCS pointing mode, which could never be activated, probably because of a failure in the star tracker's Lost in Space Mode; and another due to the electronic failure of a critical system by unknown reasons, which caused loss of communication with the CubeSat.

After ten weeks of operations, on 2018 April 5 the mission was officially ended. Although unfortunately almost no payload data besides payload temperatures were obtained, sparse but still valuable platform data were acquired Nowak et al. (2018).

A crucial step before modeling is to define the set of reference frames necessary for describing attitude and orbit kinematics (positions and velocities) as well as dynamics. We adopt the frame notation presented in Craig (1989), to avoid ambiguity about the reference frame in which a quantity is computed and the one in which it is expressed.

Consider three reference frames {A}, {B}, and {C}, with {B} rotating with velocity Ωおめが. The velocity of frame {B} with respect to frame {A} is denoted ^A Ωおめが_B , and ^C (^A Ωおめが_B ) denotes the velocity of frame {B} with respect to frame {A} expressed in frame {C}. If {A} is an inertial frame, we simplify notation by using lowercase letters. Thus, the velocity of frame {B} with respect to the inertial frame {A} is ωおめが _B , and ^C ωおめが_B is the velocity of frame {B} with respect to the inertial frame {A} expressed in frame {C}. Figure 5 illustrates the four reference frames necessary for the modeling presented in this section.

• 1.  
  Earth-Centered Inertial Reference Frame (ECI), denoted by {I}. The most important reference frame is the inertial one. All the physical quantities are computed with respect to this frame, and only afterwards expressed with respect to another useful frame. See Appendix A.1.1 for axes definition.
• 2.  
  Earth-Centered Earth-Fixed Reference Frame (ECEF), denoted by {E}, is a non-inertial frame. See Appendix A.1.2 for axes definition.
• 3.  
  Spacecraft Reference Frame (SCRF), denoted by {S}, is useful to represent some elements of the spacecraft, such as actuators, and is fixed with respect to the satellite body. See Appendix A.1.3 for axes definition.
• 4.  
  Attitude Control Reference Frame (ACRF), denoted by {A}, is useful for attitude control. The spacecraft dynamics is simplified if represented in relation to its principal axes of inertia. See Appendix A.1.4 for axes definition.

Zoom In Zoom Out Reset image size

Appendix A.2 presents the transformations between these reference frames. The next subsections present the key expressions for attitude and orbit dynamics in the reference frames. Derivations and further details are presented in Appendix.

The overall satellite attitude dynamic model is shown in Equation (1). Its derivation can be found in Appendices A.3–A.6

$$\begin{eqnarray}&&\begin{array}{l}{}^{A}{{\boldsymbol{T}}}_{\mathrm{mtq}}-{}^{A}{\dot{{\boldsymbol{H}}}}_{\mathrm{rw}}+{}^{A}{{\boldsymbol{T}}}_{\mathrm{dst}}\\ \ =\,{}^{A}{{\boldsymbol{I}}}_{\mathrm{sc}}{}^{A}{\dot{{\boldsymbol{\omega }}}}_{A}+{}^{A}{{\boldsymbol{\omega }}}_{A}\times \left({}^{A}{{\boldsymbol{I}}}_{\mathrm{sc}}{}^{A}{{\boldsymbol{\omega }}}_{A}+{}^{A}{{\boldsymbol{H}}}_{\mathrm{rw}}\right),\end{array}\end{eqnarray} \tag{ 1 }$$

where ^A T _mtq is the total magnetic actuators torque, ${}^{A}{\dot{{\boldsymbol{H}}}}_{\mathrm{rw}}$ is the controllable reaction wheels torque, ^A T _dst is the torque due to all disturbances, ^A I _sc is the spacecraft inertia matrix, ${}^{A}{\dot{{\boldsymbol{\omega }}}}_{A}$ and ^A ωおめが _A are the satellite angular acceleration and velocity, respectively. Here ^A ωおめが _A is the angular velocity of attitude satellite frame {A} with respect to inertial frame {I} and expressed in terms of frame {A}. Finally, ^A H _rw is the angular momentum of the reaction wheels.

4.2.1. Attitude Disturbance

Satellites are subject to disturbances that depend on a series of factors including mechanical characteristics of the spacecraft, mission profile and solar activity. In LEO, the main torques are due to the atmospheric (aerodynamic) drag, T _d, the solar radiation pressure, T _srp, the gradient of gravity, T _gg, and the residual magnetic dipole, T _mag. The torque due to all disturbances, T _dst, is the sum of all these torques expressed in terms of frame {S}:

$$\begin{eqnarray}&&{}^{S}{{\boldsymbol{T}}}_{\mathrm{dst}}={}^{S}{{\boldsymbol{T}}}_{{\rm{d}}}+{}^{S}{{\boldsymbol{T}}}_{\mathrm{srp}}+{}^{S}{{\boldsymbol{T}}}_{\mathrm{gg}}+{}^{S}{{\boldsymbol{T}}}_{\mathrm{mag}}.\end{eqnarray} \tag{ 2 }$$

Appendix A.7.1 details these terms.

One way to obtain orbit propagation is through Cowell's Formulation (Vallado 2013), a Keplerian orbit modified model considering a disturbance acceleration a _dst:

$$\begin{eqnarray}&&{}^{I}{{\boldsymbol{a}}}_{\mathrm{sc}}=-\displaystyle \frac{{\mu }_{\oplus }}{{{\boldsymbol{r}}}_{\mathrm{sc}}^{3}}{}^{I}{{\boldsymbol{r}}}_{\mathrm{sc}}+{}^{I}{{\boldsymbol{a}}}_{\mathrm{dst}},\end{eqnarray} \tag{ 3 }$$

where μみゅー_⊕ is the orbital parameter of the Earth equal to 398,600 km³ s⁻², r _sc is the position vector of the spacecraft relative to the Earth and a _dst is the acceleration due to all orbital disturbances written in {I}.

This modeling choice provides a high-accuracy orbit propagator as will be verified is Section 5.2.

4.3.1. Orbital Disturbances

The most important orbital disturbances in LEO are the atmospheric drag, a _d, solar radiation pressure, a _srp, Earth's oblateness, a _ob, and third body influence (Sun and Moon), a _th. The acceleration due to all disturbances, a _dst, is given by adding the disturbance terms:

$$\begin{eqnarray}&&{}^{I}{{\boldsymbol{a}}}_{\mathrm{dst}}={}^{I}{{\boldsymbol{a}}}_{{\rm{d}}}+{}^{I}{{\boldsymbol{a}}}_{\mathrm{srp}}+{}^{I}{{\boldsymbol{a}}}_{\mathrm{ob}}+{}^{I}{{\boldsymbol{a}}}_{\mathrm{th}}.\end{eqnarray} \tag{ 4 }$$

Appendix A.7.2 details these terms.

Sensor models consider measurement of noise and bias. Gyrometer measurements are expressed by

$$\begin{eqnarray}&&{\boldsymbol{\omega }}={{\boldsymbol{\omega }}}_{\mathrm{actual}}+{{\boldsymbol{b}}}_{\omega }+{{\boldsymbol{\sigma }}}_{\omega },\end{eqnarray} \tag{ 5 }$$

where ωおめが _actual is the actual angular velocity, b _ωおめが is the bias and σしぐま _ωおめが is a white noise.

Similarly, magnetometer measurements are expressed by

$$\begin{eqnarray}&&{\boldsymbol{B}}={{\boldsymbol{B}}}_{\mathrm{actual}}+{{\boldsymbol{b}}}_{{\rm{B}}}+{{\boldsymbol{\sigma }}}_{{\rm{B}}},\end{eqnarray} \tag{ 6 }$$

where B _actual is the actual Earth's magnetic field orientation, b _B is the bias and σしぐま _B is a white noise.

Since PicSat had no sensors to measure torques provoked by external disturbances, attitude model identification must be performed in four steps.

• Step 1: identification of the the open-loop satellite dynamics considering external disturbances as negligible.
• Step 2: validation of model obtained in Step 1.
• Step 3: identification of the open-loop external disturbances, given the validated satellite dynamics.
• Step 4: validation of model obtained in Step 3.

Figure 7 presents these steps in a flow diagram. One can refer to this figure during steps explanation.

Zoom In Zoom Out Reset image size

We use the satellite dynamic model presented in Section 4.2 to determine the time window regarding PicSat open-loop data for which we can consider the external disturbances negligible. Isolating the derivative of the spacecraft angular velocity in Equation (1) gives

$$\begin{eqnarray}\begin{array}{rcl}{}^{A}{\dot{{\boldsymbol{\omega }}}}_{A} & = & {\left({}^{A}{{\boldsymbol{I}}}_{\mathrm{sc}}\right)}^{-1}\left[{}^{A}{{\boldsymbol{T}}}_{\mathrm{mtq}}-{}^{A}{\dot{{\boldsymbol{H}}}}_{\mathrm{rw}}+{}^{A}{{\boldsymbol{T}}}_{\mathrm{dst}}\right.\\ & & -\left.\,{}^{A}{{\boldsymbol{\omega }}}_{A}\times \left({}^{A}{{\boldsymbol{I}}}_{\mathrm{sc}}{}^{A}{{\boldsymbol{\omega }}}_{A}+{}^{A}{{\boldsymbol{H}}}_{\mathrm{rw}}\right)\right].\end{array}\end{eqnarray} \tag{ 7 }$$

Equation (7) is the dynamic simulation form of Equation (1).

Considering the case in which magnetorquers are off ( T _mtq = 0) and reaction wheels stopped ( H _rw = 0), Equation (7) becomes

$$\begin{eqnarray}&&{}^{A}{\dot{{\boldsymbol{\omega }}}}_{A}={\left({}^{A}{{\boldsymbol{I}}}_{\mathrm{sc}}\right)}^{-1}\left[{}^{A}{{\boldsymbol{T}}}_{\mathrm{dst}}-{}^{A}{{\boldsymbol{\omega }}}_{A}\times \left({}^{A}{{\boldsymbol{I}}}_{\mathrm{sc}}{}^{A}{{\boldsymbol{\omega }}}_{A}\right)\right].\end{eqnarray} \tag{ 8 }$$

Omitting the sub and superscribed reference frames and taking the integration of Equation (8) between times t_i and t_f, we obtain

$$\begin{eqnarray}&&{\boldsymbol{\omega }}({t}_{{\rm{f}}})={\boldsymbol{\omega }}({t}_{{\rm{i}}})-{{\boldsymbol{I}}}_{\mathrm{sc}}^{-1}{\int }_{{t}_{{\rm{i}}}}^{{t}_{{\rm{f}}}}{\boldsymbol{\omega }}\times ({{\boldsymbol{I}}}_{\mathrm{sc}}{\boldsymbol{\omega }}){dt}+{{\boldsymbol{I}}}_{\mathrm{sc}}^{-1}{\int }_{{t}_{{\rm{i}}}}^{{t}_{{\rm{f}}}}{{\boldsymbol{T}}}_{\mathrm{dst}}{dt}.\end{eqnarray} \tag{ 9 }$$

The first term in Equation (9) is the angular velocity of free motion, and depends only on the satellite initial velocity and its own gyroscopic torque:

$$\begin{eqnarray}&&{{\boldsymbol{\omega }}}_{\mathrm{free}}({t}_{{\rm{f}}})={\boldsymbol{\omega }}({t}_{{\rm{i}}})-{{\boldsymbol{I}}}_{\mathrm{sc}}^{-1}{\int }_{{t}_{{\rm{i}}}}^{{t}_{{\rm{f}}}}{\boldsymbol{\omega }}\times ({{\boldsymbol{I}}}_{\mathrm{sc}}{\boldsymbol{\omega }}){dt}.\end{eqnarray} \tag{ 10 }$$

In Equation (10), if $\parallel {\boldsymbol{\omega }}({t}_{{\rm{i}}})\parallel \gg \parallel {{\boldsymbol{I}}}_{\mathrm{sc}}^{-1}{\int }_{{t}_{{\rm{i}}}}^{{t}_{{\rm{f}}}}{\boldsymbol{\omega }}\times ({{\boldsymbol{I}}}_{\mathrm{sc}}{\boldsymbol{\omega }}){dt}\parallel$ then ∥ ωおめが _free(t_f)∥ ≈ ∥ ωおめが (t_i)∥. The second term in Equation (9) is the angular velocity due to external disturbances:

$$\begin{eqnarray}&&{{\boldsymbol{\omega }}}_{\mathrm{dst}}({t}_{{\rm{f}}})={{\boldsymbol{I}}}_{\mathrm{sc}}^{-1}{\int }_{{t}_{{\rm{i}}}}^{{t}_{{\rm{f}}}}{{\boldsymbol{T}}}_{\mathrm{dst}}{dt}.\end{eqnarray} \tag{ 11 }$$

The worst case for Equation (11) occurs when T _dst is maximum and constant, that is,

$$\begin{eqnarray*}&&{{\boldsymbol{\omega }}}_{{\mathrm{dst}}_{\max }}({t}_{{\rm{f}}})={{\boldsymbol{I}}}_{\mathrm{sc}}^{-1}{{\boldsymbol{T}}}_{{\mathrm{dst}}_{\max }}({t}_{{\rm{f}}}-{t}_{{\rm{i}}}).\end{eqnarray*}$$

Thus, the ratio between the free motion and the maximum disturbance in terms of norm of angular velocity is

$$\begin{eqnarray*}&&\displaystyle \frac{\parallel {{\boldsymbol{\omega }}}_{\mathrm{free}}({t}_{{\rm{f}}})\parallel }{\parallel {{\boldsymbol{\omega }}}_{{\mathrm{dst}}_{\max }}({t}_{{\rm{f}}})\parallel }\approx \displaystyle \frac{\parallel {\boldsymbol{\omega }}({t}_{{\rm{i}}})\parallel }{\parallel {{\boldsymbol{I}}}_{\mathrm{sc}}^{-1}{{\boldsymbol{T}}}_{{\mathrm{dst}}_{\max }}\parallel ({t}_{{\rm{f}}}-{t}_{{\rm{i}}})}.\end{eqnarray*}$$

In attitude modeling for closed-loop control, we may neglect disturbances of less than one tenth of the open-loop response amplitude, that is, ∥ ωおめが _free∥/∥ ωおめが _dst∥ ≥ 10. Therefore

$$\begin{eqnarray*}&&({t}_{{\rm{f}}}-{t}_{{\rm{i}}})\leqslant \displaystyle \frac{\parallel {\boldsymbol{\omega }}({t}_{{\rm{i}}})\parallel }{10\parallel {{\boldsymbol{I}}}_{\mathrm{sc}}^{-1}{{\boldsymbol{T}}}_{{\mathrm{dst}}_{\max }}\parallel }.\end{eqnarray*}$$

Considering the suitable in-orbit data subsets for such evaluation, the lowest norm of the angular velocity was ωおめが (t_i) = 0.028 rad s⁻¹.

Taking ${{\boldsymbol{I}}}_{\mathrm{sc}}^{-1}$ from the PicSat CADきゃど model as the best initial approximation of the true ${{\boldsymbol{I}}}_{\mathrm{sc}}^{-1}$, and ${{\boldsymbol{T}}}_{{\mathrm{dst}}_{\max }}$ estimated via the attitude dynamic model using initial values for parameters taken from satellite components data sheets, we find $\parallel {{\boldsymbol{I}}}_{\mathrm{sc}}^{-1}{{\boldsymbol{T}}}_{{\mathrm{dst}}_{\max }}\parallel =4.15\times {10}^{-5}$ s⁻². Thus

$$\begin{eqnarray*}&&{\rm{\Delta }}t\leqslant \displaystyle \frac{0.028}{\mathrm{10.4.15}\times {10}^{-5}}=67.5\,{\rm{s}}.\end{eqnarray*}$$

Therefore Step 1 must be applied to time windows up to 67.5 s. Since ${{\boldsymbol{T}}}_{{\mathrm{dst}}_{\max }}$ is a very conservative estimation and a better value can be determined by iterating modeling results against in-orbit data, we will relax this constraint to 100 s allowing the use of more in-orbit data and improving model accuracy.

In order to measure ωおめが (t_i) and ${{\boldsymbol{T}}}_{{\mathrm{dst}}_{\max }}$ in open-loop attitude response with magnetorquers off and reaction wheels stopped, analysis of in-orbit angular velocity data allowed us to identify 10 data stretches with at least 10 consecutive samples having the following lengths: 290, 680, 570, 300, 770, 600, 690, 750, 720 and 620 s. To carry out Steps 1 and 2, these stretches were divided into subsets of 100 s. When necessary to complete the last subset, data from the previous one was used. There were 63 subsets altogether.

5.1.1. Step 1—Satellite Dynamics Identification

For Step 1, 32 from the original 63 subsets were selected. An arbitrary rule was applied: the odd index subsets were chosen, that is, the first, the third, the fifth, and so on. PicSat inertia matrix, ADCS gyrometers' bias and initial velocity were optimized/estimated based on a nonlinear least squares approach using those subsets, that is, they were all free parameters of the optimization algorithm. Principal inertia moments varied from about 1% in X and Y axes to 10% in Z-axis comparing CADきゃど guess and identified values. Furthermore, the gyrometer biases were in the order of 10⁻³ rad s⁻¹ and their respective values were subtracted from each axis.

In order to quantify the difference between in-orbit data and the identified attitude model, we employ the metric

$$\begin{eqnarray}&&{\int }_{{t}_{{\rm{i}}}}^{{t}_{{\rm{f}}}}{\boldsymbol{diff}}\approx \sum _{k=1}^{N}{\parallel {{\boldsymbol{\omega }}}_{\mathrm{in} \mbox{-} \mathrm{orbit}}(k)-{{\boldsymbol{\omega }}}_{\mathrm{model}}(k)\parallel }_{2},\end{eqnarray} \tag{ 12 }$$

where diff is a vector containing the ℓ²-norm of the angular velocity difference between in-orbit and modeled satellite angular velocity, N is the number of measurements acquired between t_i and t_f. As Equation (12) is a summation of modeling residues, its value increases with time for each velocity axis. More classical metrics would lead us to deal with in-orbit data noise filtering, which is an extra work without any big advantage in terms of modeling accuracy regarding closed-loop control.

Figure 8 shows the identification results for the angular velocity in terms of Equation (12). In this plot each dot is the result of summation of 100 s of ℓ²-norm data.

Zoom In Zoom Out Reset image size

The average differences are distributed around the abscissa axis, indicating residual randomness. Summation of norm is typically around 7 × 10⁻³ rad s⁻¹, which means that identification residuals are very small. We show them in more details in Step 2. Only for Z-axis, norm summation is as high as 2 × 10⁻² rad s⁻¹, which is due to the gyrometer noise.

5.1.2. Step 2—Satellite Dynamics Validation

The validation was performed with the other 31 subsets. In this case, even indexes were used, that is, the second, the fourth, the sixth and so on. An example of result for one subset is shown in Figure 9.

Zoom In Zoom Out Reset image size

One can see an excellent match between in-orbit and modeled angular velocity, except for Z-axis, which is not due to the dynamic model. Instead, it is due to high gyrometer noise. In this time frame (100 s), the satellite behaves as a rigid body in free motion. Thus, we have kinetic energy being exchanged among its three rotation axes, which means that angular velocities are coupled and are constrained by the rigid body dynamics. Thus, it is physically impossible to have Z-axis high-frequency velocity while X and Y axes move slowly. One can observe this behavior all around the subsets in Figure 8 as well. Therefore, for modeling and validation purposes, X and Y axes raw data are more reliable them Z-axis data. If we apply a low-pass filter matching the frequency domain of the satellite dynamics, they will be more consistently reliable. We have not designed and implemented such a filter, because a proper development requires detailed knowledge about the flight sensors while the raw data accuracy is already enough for dynamic modeling in terms of ADCS.

Figure 10 shows the validation results for the angular velocity in terms of Equation (12). In this plot each dot is the result of summation of 100 s of ℓ²-norm data. In other words, Figure 9 becomes three dots in Figure 10, one for each axis.

Zoom In Zoom Out Reset image size

As for the identification results, the average differences appear to be random. This time, Figure 9 helps to understand what the amplitudes mean. The summation of the norm for the X-axis is 0.5 × 10⁻² rad s⁻¹, for the Y-axis is 1 × 10⁻² rad s⁻¹ and for the Z-axis is 2 × 10⁻² rad s⁻¹. Therefore, 1 ×10⁻² rad s⁻¹ or less indicates a really good match between real data and model.

5.1.3. Step 3—External Disturbances Identification

To accommodate uncertainties on the disturbance model parameters, and even data time stamping bias, a disturbance gain, K_dst, was introduced in Equation (2), leading to:

$$\begin{eqnarray}&&{}^{S}{{\boldsymbol{T}}}_{\mathrm{dst}}={}^{S}{K}_{\mathrm{dst}}({}^{S}{{\boldsymbol{T}}}_{{\rm{d}}}+{}^{S}{{\boldsymbol{T}}}_{\mathrm{srp}}+{}^{S}{{\boldsymbol{T}}}_{\mathrm{gg}}+{}^{S}{{\boldsymbol{m}}}_{\mathrm{mag}}\times {}^{S}{\boldsymbol{B}}),\end{eqnarray} \tag{ 13 }$$

where the residual magnetic dipole ( m _mag) will be also a degree of freedom, as long as it is the most important disturbance in the specific case of PicSat.

As in Step 1, to perform Step 3 we selected 5 subsets alternately from the original 10 ones. The subsets are now longer than 100 s. The validated inertia matrix and gyro bias are obtained from Step 2. As previously, the procedure to identify K_dst and m _mag utilized a nonlinear least squares approach. The values obtained were

$$\begin{eqnarray*}&&\begin{array}{l}{}^{S}{{\boldsymbol{m}}}_{\mathrm{mag}}\\ \ =\,{\left[1.86\pm 0.04\quad 8.47\pm 0.08\quad 33.02\pm 0.65\right]}^{{\rm{T}}}\times {10}^{-3}{\rm{A}}\,{{\rm{m}}}^{2}\end{array}\end{eqnarray*}$$

and

$$\begin{eqnarray*}&&{}^{S}{{\boldsymbol{K}}}_{\mathrm{dst}}={\left[1.08\pm 0.01\quad 0.95\pm 0.01\quad 1.06\pm 0.03\right]}^{{\rm{T}}}.\end{eqnarray*}$$

Figure 11 shows the identification results for the angular velocity.

Zoom In Zoom Out Reset image size

Step 3 utilizes a small number of subsets when compared to Steps 1 and 2. Another difference is that this time the duration is as long as 700 s. Even so, the summation of norm gives results in the order of 2 to 4 × 10⁻² rad s⁻¹, which is remarkably low. Therefore, disturbances have been adequately modeled.

Feedback control can deal with disturbances which are unknown or known with some degree of accuracy. Step 3 results are sufficient for designing a compensator for PicSat.

5.1.4. Step 4—External Disturbances Validation

External disturbances validation was performed using the other 5 subsets. An example of result for one subset is shown in Figure 12.

Zoom In Zoom Out Reset image size

There is an excellent match between in-orbit and modeled angular velocity even in the presence of external disturbances. The spectral density of differences is not concentrated in any specific frequency, which indicates that no systematic modeling error is present.

Figure 13 presents a zoom-in in the modeled external disturbances. The torque due to the residual magnetic dipole, T _mag, is the most important for PicSat. It is partially a consequence of letting m _mag as a free parameter and partially because, indeed, the residual magnetic dipole is typically the most relevant external disturbance torque for CubeSats. Nevertheless, any other choice would have produced a distortion in terms of disturbance weight and would not be appropriate. Figure 13 also shows that the torque due to the solar radiation pressure, T _srp, is the less important.

Zoom In Zoom Out Reset image size

In Section 5.1, using theoretical and nominal values, we have determined that for time frames smaller than 100 s, the influence of the external disturbances would be ten times less than the free motion dynamics. To verify this premise, Figure 14 presents simulations comparing the presence and absence of disturbances in the satellite's dynamic model against in-orbit data. One can observe a clear increase in the modeling error after 150 s if disturbance is not taking into account in the model.

Zoom In Zoom Out Reset image size

More precisely, when simulation reaches 200 s, the ratio between external torque and free motion reaches about 1/10. Figure 15 shows the ratio evolution regarding the time window. Such result is computed based on the validated model.

Zoom In Zoom Out Reset image size

The time window of 100 s adopted is consistent with negligible disturbance. The whole set of validation results for angular velocity is presented in Figure 16. As in the estimation results, the summation stays at very low values indicating the modeling accuracy achieved suffices for feedback control.

Zoom In Zoom Out Reset image size

The orbital position of a satellite can be obtained by means of a GPS device on board the spacecraft, or by tracking radars on the ground. PicSat lacked either system. However, the U.S. Space Surveillance Network tracks all objects in Earth orbit (artificial satellites and debris) and maintains a cataloged database of them (Vallado 2013). Their data are published by the Joint Space Operations Center on the portal Space-Track.org ¹⁰ in a format called Two-Line Element (TLE).

To estimate its current position, the ODS has an embedded orbit propagator, ¹¹ which receives the last available TLE from the ground station.

For a given object, the time between TLE updates is not regular. It can vary from a fraction of the orbit to several days. Figure 17 shows the histogram of the update time between the PicSat TLEs from 2018 January 12 through 2019 September 29.

Zoom In Zoom Out Reset image size

Altogether there were 1166 updates during this time span. The longest interval between TLEs was 30 orbits, equivalent to 47.4 hr and 75% of updates occurred in less than a 11 hr interval.

The published TLEs lack accuracy information. However, in Riesing (2015) the Planet CubeSats constellation Flock 1B was used to determine the accuracy of the TLEs of these spacecraft, once they were equipped with a GPS module. As with PicSat, they were 3Us CubeSats, launched from the ISS. Ten satellites were deployed during 2014 September and 634 TLE updates occurred during that time. The statistical error in position could be estimated with good statistical accuracy to be less than 2.01 km in 25% of cases (first quartile), less than 4.52 km in 50% of cases (median), and less than 10.60 km in 75% of cases (third quartile).

Although contact with PicSat was lost on 2018 March 20, TLEs will continue to be available until the satellite re-enters Earth's atmosphere. We performed the identification and validation of the simulator's orbital propagation using the available TLEs for PicSat and the statistical error obtained in Riesing (2015). For this purpose, the following procedure was carried out: using TLEs as initial conditions for the simulator, we propagate the states until the next available TLE, comparing them at the end. We repeat this process ranging from TLEs available between one orbit up to 16 orbits (i.e., about 1 day between two TLEs). The results of the propagation error are shown in Figure 18. The horizontal dashed lines represent the statistical errors obtained in Riesing (2015).

Zoom In Zoom Out Reset image size

Model error increases with the propagation time. After 15 orbits, the median error of the model is 7.1 km, which is a very small error considering that 648,047.7 km were covered, i.e., 1.1 × 10⁻⁷% propagation error in 15 orbits. For an observer on the ground, it is equivalent to $3\buildrel{\,\prime}\over{.} 5$ uncertainty in pointing toward the spacecraft. For payloads that do not demand very high accuracy orbit determination, such propagator accuracy is sufficient.

Based on this model, the predicted re-enter date for PicSat is 2022 December 24 03:44:34 (UTC).

Figure 19 shows the PicSat angular velocity norm during the entire mission. Start and stop commands for de-tumbling mode are signaled with dashed and dotted vertical lines, respectively. Three of the occasions, where an effective reduction in angular velocity can be seen, are marked by numbers 1–3.

Zoom In Zoom Out Reset image size

Considering the B-dot control structure, ${{\boldsymbol{k}}}_{\det }$ gain and initial velocity were estimated in a way similar to the dynamics validation performed in Section 5.1, using Equation (12). Free parameters of the algorithm were the initial velocity and the controller gain. Once again, the procedure was performed using a nonlinear least squares method. Data stretch number 2 was used to do the estimation and data stretches number 1 and 3 were used for the validation. B-dot gain identification and validation results are shown in Figure 20.

Zoom In Zoom Out Reset image size

The identified controller (Equation (14)) presents a satisfactory result for the norm of angular velocities. This indicates that the kinetic energy reduction is achieved in a similar way for both flight and on-ground modeled controllers. However, if the in-orbit and modeled angular velocities are compared independently, axis by axis, the results are unsatisfactory, suggesting some inconsistency between them. Such an issue will be investigated and an advanced identification strategy as discussed in Romano & Pait (2017) will be applied to precisely determine what is missing in Equation (14).

We start by clearly defining all necessary reference frames.

A.1.1. Earth-centered Inertial Reference Frame (ECI)

A.1.2. Earth-centered Earth-fixed Reference Frame (ECEF)

A.1.3. Spacecraft Reference Frame (SCRF)

A.1.4. Attitude Control Reference Frame (ACRF)

Reference frames can be related to each other by frame transformations.

A.2.1. Transformation between ECI and ECEF

Once the ECI and ECEF ($\hat{X}\hat{Y}$)-planes are coplanar, the rotation between them comes only from the Earth's rotation movement around the $\hat{Z}$ axis. Thus, the rotation matrix and its representation in quaternions are expressed by

$$\begin{eqnarray}{}^{I}{{\boldsymbol{Q}}}_{E}=\left[\begin{array}{ccc}\cos {\theta }_{{\rm{G}}} & \sin {\theta }_{{\rm{G}}} & 0\\ -\sin {\theta }_{{\rm{G}}} & \cos {\theta }_{{\rm{G}}} & 0\\ 0 & 0 & 1\end{array}\right]\end{eqnarray} \tag{ A1 }$$

and

$$\begin{eqnarray}&&{}^{I}{q}_{E}={\left[\cos \left(\displaystyle \frac{{\theta }_{{\rm{G}}}}{2}\right)\quad 0\quad 0\quad \sin \left(\displaystyle \frac{{\theta }_{{\rm{G}}}}{2}\right)\right]}^{{\rm{T}}},\end{eqnarray} \tag{ A2 }$$

where θしーた_G is the Greenwich Sidereal Time, that is, the angle between the March equinox and the main meridian in degrees. As shown in Vallado (2013),

$$\begin{eqnarray}&&{\theta }_{{\rm{G}}}={\theta }_{{\rm{G}}0}+360.98564724\,\displaystyle \frac{{t}_{{\rm{U}}}}{24},\end{eqnarray} \tag{ A3 }$$

where t_U is the time in hours in UT1 scale, and θしーた_G0 is the Greenwich Sidereal Time at 0 hr UT1. As presented in Seidelmann (2006), θしーた_G0 is given by

$$\begin{eqnarray}\begin{array}{rcl}{\theta }_{{\rm{G}}0} & = & 100.4606184+36000.77004\,{T}_{0}\\ & & +\,0,000387933\,{T}_{0}^{2}-2.583\times {10}^{-8}\,{T}_{0}^{3},\end{array}\end{eqnarray} \tag{ A4 }$$

where T₀ = d_U/36525, and d_U is the time in Julian days since J2000.0.

A.2.2. Transformation between SCRF and ACRF

The transformation of the spacecraft inertia matrix I _sc, from SCRF or just {S} to ACRF or just {A}, is given by

$$\begin{eqnarray}&&{}^{A}{\boldsymbol{I}}=\left({}^{A}{{\boldsymbol{Q}}}_{S}\right)\left({}^{S}{\boldsymbol{I}}\right){\left({}^{A}{{\boldsymbol{Q}}}_{S}\right)}^{{\rm{T}}},\end{eqnarray} \tag{ A5 }$$

where ^A I is the diagonal form of the spacecraft's inertia matrix, ^S I is the non-diagonal form of the spacecraft's inertia matrix and ^A Q _S is the rotation matrix from {S} to {A}.

As shown in Curtis (2014), the matrix ^A I is formed, on the main diagonal, by the eigenvalues of ^S I and the matrix ^A Q _S is formed by the respective eigenvectors of ^S I .

We can represent the satellite orientation by means of rotation matrices (also called direction cosine matrices), Euler angles, or quaternions. Rotation matrices are hard to read and computationally expensive. Euler angles are easy to read but present singularities and require the use of rotation matrices when transforming reference frames. Quaternions are hard to read and computationally inexpensive. Therefore, from now on we present results using Euler angles, and compute transformations using quaternions, which are preferable for embedded ACS (Sidi 1997; Curtis 2014).

The Euler angles are three angles describing the orientation of a rigid body with respect to a fixed coordinate system. Here they are expressed by (ϕ, θしーた, ψぷさい) corresponding to yaw, pitch and roll angles in frame {A}.

A quaternion can be expressed as:

$$\begin{eqnarray}\begin{array}{rcl}q & = & {q}_{0}+{q}_{1}\hat{{\boldsymbol{i}}}+{q}_{2}\hat{{\boldsymbol{j}}}+{q}_{3}\hat{{\boldsymbol{k}}}\\ & = & {q}_{0}+{\boldsymbol{q}};\end{array}\end{eqnarray} \tag{ A6 }$$

where q₀ is the scalar part and q is the vector or imaginary part (Wertz 1978).

By restricting the quaternion to be unitary, $\parallel q\parallel =\sqrt{{q}_{0}^{2}+{q}_{1}^{2}+{q}_{2}^{2}+{q}_{3}^{2}}=1$, it is possible to relate it to the symmetrical Euler parameters.

The transformation between them using the 3-2-1 sequence is given by:

$$\begin{eqnarray}&&{}^{I}{q}_{A}=\left[\begin{array}{c}\cos \left(\displaystyle \frac{\psi }{2}\right)\cos \left(\displaystyle \frac{\theta }{2}\right)\cos \left(\displaystyle \frac{\phi }{2}\right)+\sin \left(\displaystyle \frac{\psi }{2}\right)\sin \left(\displaystyle \frac{\theta }{2}\right)\sin \left(\displaystyle \frac{\phi }{2}\right)\\ \sin \left(\displaystyle \frac{\psi }{2}\right)\cos \left(\displaystyle \frac{\theta }{2}\right)\cos \left(\displaystyle \frac{\phi }{2}\right)-\cos \left(\displaystyle \frac{\psi }{2}\right)\sin \left(\displaystyle \frac{\theta }{2}\right)\sin \left(\displaystyle \frac{\phi }{2}\right)\\ \cos \left(\displaystyle \frac{\psi }{2}\right)\sin \left(\displaystyle \frac{\theta }{2}\right)\cos \left(\displaystyle \frac{\phi }{2}\right)+\sin \left(\displaystyle \frac{\psi }{2}\right)\cos \left(\displaystyle \frac{\theta }{2}\right)\sin \left(\displaystyle \frac{\phi }{2}\right)\\ \cos \left(\displaystyle \frac{\psi }{2}\right)\cos \left(\displaystyle \frac{\theta }{2}\right)\sin \left(\displaystyle \frac{\phi }{2}\right)-\sin \left(\displaystyle \frac{\psi }{2}\right)\sin \left(\displaystyle \frac{\theta }{2}\right)\cos \left(\displaystyle \frac{\phi }{2}\right)\end{array}\right].\end{eqnarray} \tag{ A7 }$$

Concerning the satellite, kinematics describes its attitude and angular velocity time series without any concern about the torques that have created such movement.

Consider the spacecraft's Cartesian angular velocity described by

$$\begin{eqnarray}&&{{\boldsymbol{\omega }}}_{A}={\left[{\omega }_{x}\quad {\omega }_{y}\quad {\omega }_{z}\right]}_{A}^{{\rm{T}}},\end{eqnarray} \tag{ A8 }$$

where ωおめが_x , ωおめが_y and ωおめが_z are, respectively, the angular velocities around x, y, and z axes of frame {A} w.r.t. frame {I}. As shown in Sidi (1997), the relationship between the time derivative of quaternions and the Cartesian angular velocity is given by:

$$\begin{eqnarray}&&\displaystyle \frac{d}{{dt}}{}^{I}{q}_{A}=\displaystyle \frac{1}{2}{\boldsymbol{W}}({{\boldsymbol{\omega }}}_{A}){}^{I}{q}_{A},\end{eqnarray} \tag{ A9 }$$

where W ( ωおめが _A ) is a skew-symmetric matrix given by

$$\begin{eqnarray}{\boldsymbol{W}}({{\boldsymbol{\omega }}}_{A})={\left[\begin{array}{cccc}0 & {\omega }_{z} & -{\omega }_{y} & {\omega }_{x}\\ -{\omega }_{z} & 0 & {\omega }_{x} & {\omega }_{y}\\ {\omega }_{y} & -{\omega }_{x} & 0 & {\omega }_{z}\\ -{\omega }_{x} & -{\omega }_{y} & -{\omega }_{z} & 0\end{array}\right]}_{A}.\end{eqnarray} \tag{ A10 }$$

As ωおめが_x , ωおめが_y and ωおめが_z are measured with respect to the inertial frame {I}, the satellite attitude is simply their integral. Knowing the satellite's angular velocity, it is straightforward to compute the equivalent quaternion derivative and the quaternion itself, which is useful for simulation purposes.

Attitude dynamics describes the angular velocity as a function of the applied torque T .

Most CubeSats can be considered rigid bodies, because they do not have moving parts and orbit control systems, which usually employ of propellants. Therefore, mass distribution remains constant over time. Although PicSat has a piezoelectric positioning system as part of the payload (Piezo-stage in Figure 4), it can be considered a rigid body for attitude dynamics modeling. As a result, rigid body dynamics properties and equations apply. Considering ωおめが _A the angular velocity of {ACRF} or just {A} with respect to {ECI} or just {I}, the Euler equation is given by Curtis (2014); Sidi (1997):

$$\begin{eqnarray}&&{}^{A}{{\boldsymbol{T}}}_{\mathrm{act}}+{}^{A}{{\boldsymbol{T}}}_{\mathrm{dst}}={}^{A}{{\boldsymbol{I}}}_{\mathrm{sc}}{}^{A}{\dot{{\boldsymbol{\omega }}}}_{A}+{}^{A}{{\boldsymbol{\omega }}}_{A}\times \left({}^{A}{{\boldsymbol{I}}}_{\mathrm{sc}}{}^{A}{{\boldsymbol{\omega }}}_{A}\right)\end{eqnarray} \tag{ A11 }$$

where ^A I _sc and ^A ωおめが _A are the spacecraft inertia matrix and the angular velocity written in {A}, respectively, ^A T _act is the torque provoked by actuators and ^A T _dst is the torque provoked by disturbances.

According to Table 1, PicSat is equipped with magnetic actuators and reaction wheels. Therefore, the torque due to actuators is given by Equation (A12)

$$\begin{eqnarray}&&{}^{A}{{\boldsymbol{T}}}_{\mathrm{act}}={}^{A}{{\boldsymbol{T}}}_{\mathrm{mtq}}+{}^{A}{{\boldsymbol{T}}}_{\mathrm{rw}}\end{eqnarray} \tag{ A12 }$$

where ^A T _mtq is the total torque due to magnetorquers and ^A T _rw is the total torque due to reaction wheels written in {A}.

We compute each of these torques in the next subsections.

A.6.1. Magnetorquers

In Low Earth Orbit (LEO), Earth's magnetic field intensity is not negligible (Wertz 1978). Its interaction with magnetorquers is capable of producing torque in the spacecraft. Physically, the actuator consists of a coil with a ferromagnetic core or air core. Therefore, it is possible, by Kirchhoff's Laws, to model it as a resistor and inductor series circuit. The magnetic moment is in the direction of the normal to the coil plane and the sense is given by the Right-Hand Rule. Its intensity is given by the Ampere's law. We neglect the demagnetization effect. The total intensity of the magnetic moment, m_mtq, is the sum of the intensities of the magnetic moment generated by the coil and by the ferromagnetic core. In vector form, the magnetorquer's torque is given by Equation (A13)

$$\begin{eqnarray}&&{}^{A}{{\boldsymbol{T}}}_{\mathrm{mtq}}=\sum _{i=1}^{3}{}^{A}{{\boldsymbol{m}}}_{{\mathrm{mtq}}_{i}}\times {}^{A}{\boldsymbol{B}}={}^{A}{{\boldsymbol{m}}}_{\mathrm{mtq}}\times {}^{A}{\boldsymbol{B}}\end{eqnarray} \tag{ A13 }$$

where ${{\boldsymbol{m}}}_{{\mathrm{mtq}}_{i}}$ is the magnetic moment of each dipole, and B the Earth's magnetic field.

A.6.2. Reaction Wheels

Reaction wheels (RW) allow more precise pointing control (Yost 2020). Considering I _rw and ωおめが _rw, the moment of inertia and the angular velocity of reaction wheels, and applying the Euler characteristic equation, we obtain, for each reaction wheel i:

$$\begin{eqnarray}&&{}^{A}{{\boldsymbol{T}}}_{{\mathrm{rw}}_{i}}={}^{A}{\dot{{\boldsymbol{H}}}}_{{\mathrm{rw}}_{i}}+{}^{A}{{\boldsymbol{\omega }}}_{A}\times {}^{A}{{\boldsymbol{H}}}_{{\mathrm{rw}}_{i}}.\end{eqnarray} \tag{ A14 }$$

Only the first term of Equation (A14) (${}^{A}{\dot{{\boldsymbol{H}}}}_{{\mathrm{rw}}_{i}}$) can be used to control the spacecraft orientation. The second one (${}^{A}{{\boldsymbol{\omega }}}_{A}\times {}^{A}{{\boldsymbol{H}}}_{{\mathrm{rw}}_{i}}$) is a gyroscopic perturbation effect and it will be added to the satellite dynamics. Moreover, by action and reaction law, both terms must be negative. The final equation is:

$$\begin{eqnarray}&&\begin{array}{l}{}^{A}{{\boldsymbol{T}}}_{\mathrm{mtq}}-{}^{A}{\dot{{\boldsymbol{H}}}}_{\mathrm{rw}}+{}^{A}{{\boldsymbol{T}}}_{\mathrm{dst}}={}^{A}{{\boldsymbol{I}}}_{\mathrm{sc}}{}^{A}{\dot{{\boldsymbol{\omega }}}}_{A}\\ \quad +\,{}^{A}{{\boldsymbol{\omega }}}_{A}\times \left({}^{A}{{\boldsymbol{I}}}_{\mathrm{sc}}{}^{A}{{\boldsymbol{\omega }}}_{A}+{}^{A}{{\boldsymbol{H}}}_{\mathrm{rw}}\right).\end{array}\end{eqnarray} \tag{ A15 }$$

The simplicity of Equation (A15) is suitable for PicSat modeling purposes. Nonlinear effects such as RW's imbalance (microvibrations), Coulomb friction and magnetic modeling, magnetorquer's demagnetization could be included in additional terms, depending on the requirements for each mission.

This appendix section presents all relevant attitude and orbit external disturbances that PicSat faces during flight. Figure 21 presents a block diagram that summarizes them.

Zoom In Zoom Out Reset image size

A.7.1. Attitude Disturbances

References Wertz (1978), Vallado (2013) and Curtis (2014) were used to implement attitude disturbances dynamic equations.

Torque Due to Atmospheric Drag—In orbit, the surface of a spacecraft collides with atmospheric particles, creating a force that will produce a torque in relation to the satellite's center of mass. The aerodynamic force on a surface element dA, written with respect to {S}, is given by

$$\begin{eqnarray}&&{}^{S}d{{\boldsymbol{F}}}_{{\rm{d}}}=-\displaystyle \frac{1}{2}{C}_{{\rm{D}}}{\rho }_{\mathrm{rel}}^{2}\left({}^{S}{\hat{{\boldsymbol{n}}}}_{{dA}}\cdot {}^{S}{\hat{{\boldsymbol{v}}}}_{\mathrm{rel}}\right){}^{S}{\hat{{\boldsymbol{v}}}}_{\mathrm{rel}}{dA},\end{eqnarray} \tag{ A16 }$$

where C_D is the drag coefficient, ρろー is the atmospheric density, v_rel is the spacecraft's speed relative to the atmosphere, whose direction is represented by ${\hat{{\boldsymbol{v}}}}_{\mathrm{rel}}$ and ${\hat{{\boldsymbol{n}}}}_{{dA}}$ is the unit vector of the surface element.

The torque due to atmospheric drag is expressed in {S} as

$$\begin{eqnarray}&&{}^{S}{{\boldsymbol{T}}}_{{\rm{d}}}=\int \left({}^{S}{{\boldsymbol{r}}}_{{dA}}\times {}^{S}d{{\boldsymbol{F}}}_{{\rm{d}}}\right){dA},\end{eqnarray} \tag{ A17 }$$

where r _dA is the position vector of the spacecraft's center of mass up to the surface element dA where the force d F _d is being applied.

Atmospheric Model—Density (ρろー) and other atmospheric properties vary depending on altitude, temperature, molecular composition, solar activity, etc. One of the most modern models used by the scientific community is the United States Naval Research Laboratory Mass Spectrometer and Incoherent Scatter Radar Exosphere (NRLMSISE-00) (Picone et al. 2002). It is an empirical model of the class Mass Spectrometer Incoherent Scatter Radar released in 2001 and valid up to an altitude of 1000 km. In this paper, the NRLMSISE-00 model integrated with the Simulink software is used.

Torque Due to Solar Radiation Pressure—The incidence of electromagnetic radiation on the spacecraft's surface generates disturbance torque. The main source of radiation is the direct radiation from the Sun. The intensity and direction of incidence, the shape and the optical properties of the satellite's surface are the main factors that influence the torque due to the solar radiation pressure.

Analogously to the torque due to atmospheric drag, we divide the spacecraft's surface into infinitesimal element and determine the force applied to each of them. Depending on the optical properties of the surface, radiation can be absorbed, reflected with the same incidence slope or reflected in a diffuse way. The differential radiation force due to each of these portions is given, respectively, by

$$\begin{eqnarray}&&{}^{S}d{{\boldsymbol{F}}}_{\mathrm{ra}}=\displaystyle \frac{S}{c}\left[{c}_{\mathrm{ra}}\left(-\cos \theta \ {}^{S}{\hat{{\boldsymbol{n}}}}_{{dA}}+\sin \theta \ {}^{S}{\hat{{\boldsymbol{s}}}}_{{dA}}\right)\right]\cos \theta \,{dA},\end{eqnarray} \tag{ A18 }$$

$$\begin{eqnarray}\begin{array}{rcl}{}^{S}d{{\boldsymbol{F}}}_{\mathrm{rr}} & = & \displaystyle \frac{S}{c}\left[\left(-1+{c}_{\mathrm{rr}}\right)\cos \theta \ {}^{S}{\hat{{\boldsymbol{n}}}}_{{dA}}\right.\\ & & +\,\left.\left(1-{c}_{\mathrm{rr}}\right)\sin \theta \ {}^{S}{\hat{{\boldsymbol{s}}}}_{{dA}}\right]\cos \theta \,{dA},\end{array}\end{eqnarray} \tag{ A19 }$$

and

$$\begin{eqnarray}&&{}^{S}d{{\boldsymbol{F}}}_{\mathrm{rd}}=\displaystyle \frac{S}{c}\left[-\left(\cos \theta +\displaystyle \frac{2}{3}{c}_{\mathrm{rd}}\right){}^{S}{{\boldsymbol{n}}}_{{dA}}+\sin \theta \,{}^{S}{{\boldsymbol{s}}}_{{dA}}\right]\cos \theta \,{dA},\end{eqnarray} \tag{ A20 }$$

where S is the radiation intensity, c is the speed of light, c_ra, c_rr and c_rd are, respectively, the absorbed, reflected and diffuse radiation coefficients, constrained to 0 ≤ c_ra + c_rr + c_rd ≤ 1. θしーた is the angle of incidence in relation to the normal, ${\hat{{\boldsymbol{n}}}}_{{dA}}$ and ${\hat{{\boldsymbol{s}}}}_{{dA}}$ are, respectively, the unit vector normal to the surface element and the unit vector parallel to the surface element and in the direction of radiation projection.

Thus the torque due to solar radiation pressure is given by

$$\begin{eqnarray}\begin{array}{rcl}{}^{S}{{\boldsymbol{T}}}_{\mathrm{rad}} & = & \int \left[\left({}^{S}{{\boldsymbol{r}}}_{{dA}}\times {}^{S}d{{\boldsymbol{F}}}_{\mathrm{ra}}\right)+\left({}^{S}{{\boldsymbol{r}}}_{{dA}}\times {}^{S}d{{\boldsymbol{F}}}_{\mathrm{rr}}\right)\right.\\ & & +\,\left.\left({}^{S}{{\boldsymbol{r}}}_{{dA}}\times {}^{S}d{{\boldsymbol{F}}}_{\mathrm{rd}}\right)\right]{dA},\end{array}\end{eqnarray} \tag{ A21 }$$

where r _dA is the position vector of the spacecraft's center of mass up to the surface element dA where forces d F are being applied.

Torque Due to Gravity Gradient—The gravitational force depends on the distance to the object and its orientation. Therefore, each mass element of the spacecraft is subject to a different force. This force distribution produces torque when the mass of the satellite is not perfectly distributed. This phenomenon is referred to as the gravity gradient.

The torque due to the gravity gradient, ^A T _gg, when is written with respect to {A}, is given by

$$\begin{eqnarray}&&{}^{A}{{\boldsymbol{T}}}_{\mathrm{gg}}=\displaystyle \frac{3{\mu }_{\oplus }}{{r}_{\mathrm{sc}}^{3}}\left[{}^{A}{\hat{{\boldsymbol{r}}}}_{\mathrm{sc}}\times \left({}^{A}{{\boldsymbol{I}}}_{\mathrm{sc}}{}^{A}{\hat{{\boldsymbol{r}}}}_{\mathrm{sc}}\right)\right]\end{eqnarray} \tag{ A22 }$$

where μみゅー_⊕ is the orbital parameter of the Earth, r _sc is the position vector of the spacecraft relative to the Earth and I _sc is the spacecfraft inertia matrix.

Torque Due to Residual Magnetic Dipole—The interaction between the Earth's magnetic field and the spacecraft's magnetic dipole produces torque. This magnetic dipole may be desired, like the magnetic actuator presented in Appendix A.6.1, or undesired, like eddy current, hysteresis and other magnetic moments that may arise due to currents circulating in the electronic circuits of the spacecraft. Calling it the residual magnetic dipole, m _mag, and using Equation (A13), the torque due to the residual magnetic dipole, T _mag, is given by

$$\begin{eqnarray}&&{}^{A}{{\boldsymbol{T}}}_{\mathrm{mag}}={}^{A}{{\boldsymbol{m}}}_{\mathrm{mag}}\times {}^{A}{\boldsymbol{B}},\end{eqnarray} \tag{ A23 }$$

in which B is the Earth's magnetic field and its modeling is treated below.

Model of the Earth's Magnetic Field—With the help of satellites and observatories, mathematical models of the Earth's magnetic field have been developed since the 1960s. The models mostly used by the scientific community are the International Geomagnetic Reference Field (IGRF) (Thébault et al. 2015), developed by the International Association of Geomagnetism and Aeronomy, and the World Magnetic Model (WMM) (Chulliat et al. 2015) developed by the United State National Geophysical Data Center and British Geological Survey. Both models are updated every five years and are publicly available. Software such as MATLAB and Simulink have these models integrated. This paper uses the IGRF-12 version from 2015.

A.7.2. Orbital Disturbances

References Vallado (2013), Curtis (2014) and Markley & Crassidis (2014) were used to implement orbital disturbances dynamic equations.

Acceleration Due to Atmospheric Drag—Atmospheric drag was initially explored in Appendix A.7.1. It interferes with the spacecraft's orbit by changing its acceleration. This disturbance in acceleration, a _d, is given by

$$\begin{eqnarray}&&{}^{I}{{\boldsymbol{a}}}_{{\rm{d}}}=\displaystyle \frac{{}^{I}{{\boldsymbol{F}}}_{{\rm{d}}}}{{m}_{\mathrm{sc}}},\end{eqnarray} \tag{ A24 }$$

where m_sc is the mass of the spacecraft and F _d is the drag force as shown in Equation (A16).

Acceleration Due to Solar Radiation Pressure—Analogously to the drag, the solar radiation pressure interferes in the attitude and also in the orbit of the satellite.

The disturbance a _srp is given by

$$\begin{eqnarray}&&{}^{I}{{\boldsymbol{a}}}_{\mathrm{srp}}=\displaystyle \frac{{}^{I}{{\boldsymbol{F}}}_{\mathrm{ra}}+{}^{I}{{\boldsymbol{F}}}_{\mathrm{rr}}+{}^{I}{{\boldsymbol{F}}}_{\mathrm{rd}}}{{m}_{\mathrm{sc}}},\end{eqnarray} \tag{ A25 }$$

where m_sc is the mass of the spacecraft, F _ra is the force due to the absorbed radiation (Equation (A18)), F _rr is the force due to reflected radiation (Equation (A19)), and F _rd is the force due to diffuse radiation (Equation (A20)).

Acceleration Due to Earth's Oblateness—The Earth's gravitational potential can be obtained through spherical harmonic expansion, generating an infinite series with associated Legendre polynomials.

To calculate the Earth's oblateness acceleration a _ob, it is necessary to truncate the series to a certain degree m and order n, both parameters from Legandre polynomials. The most common implementation is performed with m = 0 and n = 2. However, as presented by Vallado (2013), for a high-fidelity mathematical model, it is necessary to consider truncation at least at m = n = 24.

The most complete gravitational model is the Earth Gravitational Model 2008 (EGM2008) published by the U.S. National Geospatial-Intelligence Agency which allows truncation up to degree 2190 and order 2159 (Pavlis et al. 2012). In this paper, EGM2008 was used and the simulations were performed with truncation at m = n = 70.

Acceleration Due to Third-Body—Kepler's Orbit takes into account the gravitational interaction between two bodies. Such a situation is often referred to as the Two-Body Problem. However, other celestial bodies, such as the Moon and the Sun, will also interact gravitationally with the spacecraft. This interaction can be considered as a disturbance and modeled as a three-body system. Applying Kepler's law and considering the Moon and the Sun as third-bodies, the acceleration is given by

$$\begin{eqnarray}&&{}^{I}{{\boldsymbol{a}}}_{\mathrm{th}}={\mu }_{{\rm{m}}}\left(\displaystyle \frac{{}^{I}{{\boldsymbol{r}}}_{{\rm{m}}/\mathrm{sc}}}{{r}_{{\rm{m}}/\mathrm{sc}}^{3}}-\displaystyle \frac{{}^{I}{{\boldsymbol{r}}}_{{\rm{m}}}}{{r}_{{\rm{m}}}^{3}}\right)+{\mu }_{\odot }\left(\displaystyle \frac{{}^{I}{{\boldsymbol{r}}}_{\odot /\mathrm{sc}}}{{r}_{\odot /\mathrm{sc}}^{3}}-\displaystyle \frac{{}^{I}{{\boldsymbol{r}}}_{\odot }}{{r}_{\odot }^{3}}\right),\end{eqnarray} \tag{ A26 }$$

where μみゅー_m and μみゅー_⊙ are, respectively, the orbital parameters of the Moon and the Sun, r _m/sc and r _⊙/sc are, respectively, the position vectors of the Moon and the Sun with respect to the spacecraft, r_m and r_⊙ are, respectively, the position vectors of the Moon and the Sun with respect to the Earth.