Abstract
A detumbling algorithm is developed to yield three-axis magnetic stabilization of a CubeSat deployed with unknown RAAN, orbit phase angle, inclination, attitude, and angular rate. Data from a three-axis magnetometer are the only input to determine both the control torque and the angular rate of the spacecraft. The algorithm is designed to produce a magnetic dipole moment which is constantly orthogonal to the geomagnetic field vector, independently of both the attitude and the angular rate of the rigid spacecraft. The angular rates are calculated in real time from magnetometer data, and the use of a second-order low-pass filter allows to rapidly reduce the measurement error within ±0.2 deg/sec. Numerical validation of the algorithm is performed, and a variety of feasible scenarios is simulated assuming the CubeSat to operate in low Earth orbit. The robustness of the algorithm, with respect to unknown deployment conditions, different sampling rates, and uncertainties on the moments of inertia of the CubeSat, is verified.
1. Introduction
The stabilization of a spacecraft after deployment is never a trivial issue, and several satellites have been lost due to anomalies or dysfunctions occurred during this phase. Typically, the angular rates at the deployment are much higher than those desired for attitude maneuvering and the satellite is said to be tumbling. Therefore, a specific detumbling control must be designed to stabilize the spacecraft within the minimum time compatible with the mission requirements. After the spin motion of the spacecraft has been damped to the desired level, the control policy can be switched to fine pointing or attitude maneuvering and mission operations can start.
Additional issues often arise when implementing such a detumbling control on a CubeSat, mainly because of the limited capabilities, when compared to bigger satellites, of the Attitude Determination and Control Systems (ADCS) commonly used. Nevertheless, CubeSat technology has grown dramatically fast over the last ten years and is nowadays opening to scientific missions. Therefore, improving their reliability, and in particular by increasing the efficiency of their ADCS, is for sure a main target for the years to come.
In this work, we propose a solution for the contingent scenario in which the ADCS sensor addressed to angular rate measurement (i.e., a rate gyro) is not capable of providing any suitable information. This can happen after a failure of the sensor itself or if the CubeSat angular rates exceed the measuring range of the sensor, which is therefore not suitable to produce accurate measurements.
Focusing on CubeSats operating in low Earth orbit (LEO), an algorithm producing both the detumbling and angular rate determination using only three orthogonal magnetorquers and the measurements by a three-axis magnetometer is developed. It is worth to outline that these two results are independent one from each other. The algorithm is specifically designed to be implemented on a low computationally efficient CubeSat onboard computer. Additional constraints are represented by the limited energy budget, the strict limit on the peak current, and the low sampling rate.
The detumbling module of the algorithm is based on the popular B-dot control, here rearranged to explicitly generate a magnetic dipole moment constantly orthogonal to the geomagnetic field vector. The stability of the control was verified numerically, simulating the detumbling of the 3 U CubeSat Tigrisat. Numerical simulations also proved that an accurate selection of the proportional control parameter results in a commanded magnetic dipole moment which never exceeds a maximum value here fixed to 0.3 Am2.
The problem of magnetometer-only attitude determination has been extensively examined through the last thirty years, with relevant results obtained first by Natanson et al., for spacecraft rotating with constant angular rate and later extended to the case with no a priori knowledge of the spacecraft state with uses of a Kalman filter [1, 2]. As reviewed by Hajiyev and Guler [3], several algorithms include a single- or multiple-step extended Kalman filter [4–9] or an unscented Kalman filter [10]. Other approaches are based on deterministic two vector methods [11–13] and are typically preferred when computational efficiency of the onboard computer is limited and lower measurement accuracy is acceptable.
In the proposed algorithm, measurements of the angular rates are based on the geometric properties relating three consecutive samples of the geomagnetic field vector, processed in real time through a second-order low-pass filter. Simulations proved that the algorithm is suitable for measuring the angular rates within a steady state error of ±0.2 deg/sec and robust with respect to unknown deployment conditions and uncertainties on the inertial properties of the CubeSat.
In Sections 3 and 4, the design for the detumbling control and the angular rate determination algorithm are presented and numerically validated, simulating the detumbling of the 3 U CubeSat Tigrisat. In Section 5, a variety of detumbling scenarios, including unknown deployment parameters and uncertainties, is simulated to evaluate the performance of the algorithm and the results are discussed.
2. Mathematical Model
Some preliminary considerations are worth being outlined to clearly define the framework in which the detumbling problem is here studied. The algorithm is developed for a CubeSat deployed in LEO, without a priori knowledge about its RAAN (), orbit phase angle (), inclination (i), attitude, and angular rate. The ADCS only includes one three-axis magnetometer and three magnetorquers mutually orthogonal and aligned with the principal axes of inertia, which define the body fixed reference frame ().
Attitude dynamics is modeled assuming the spacecraft as a rigid body on which only the magnetic dipole torque acts. This is the torque due to the interaction of the magnetic dipole moment produced by magnetorquers with the geomagnetic field vector. The nonlinear angular rate dynamics can be expressed by the following equation [14]: where is the angular rate vector, is the tensor of inertia, is the magnetic dipole moment, and is the geomagnetic field vector in . Both and its rate represent the angular velocity and acceleration of with respect to , namely, the Geo Centric Inertial (GCI) reference frame [15]. The time derivatives of the geomagnetic field vector in the two frames are related by the following equation:
Numerical simulations, whose results are discussed in Sections 3, 4, and 5, were run integrating the full nonlinear spacecraft attitude dynamics equations [16], using the fixed step solver ODE8. The values of during the time of the simulation are obtained from the International Geomagnetic Reference Field (IGRF) model.
3. Detumbling Control
The proposed detumbling control represents a variation of the classical B-dot, in which the variable is defined to be constantly orthogonal to [17]. Such a result is obtained assuming that is negligible with respect to the other variations. This approximation is reasonably accurate only if the angular rates are higher than the maximum rate of change of the geomagnetic field vector, equal to twice the orbital rate () of the spacecraft [14]. The scenarios investigated in this work verify the mentioned condition (see Section 5), and thus (2) can be rearranged as follows: where is the projection of in the direction orthogonal to .
The cross-product does not allow the determination of by inverting (3). Consequently, a guessed expression for is considered as follows:
The suitability of (4) follows from the following:
Recalling that in LEO , the rhs of (5) is approximately equal to , as it is required by (3).
The magnetic dipole moment is now expressed as follows: where is a constant parameter, namely, the control gain, to be defined in accord to the specifics of the mission. Equation (6) clearly shows that is constantly orthogonal to .
The control algorithm was validated through numerical analysis in MATLAB. Simulation parameters are those of the 3 U CubeSat Tigrisat, launched by the School of Aerospace Engineering in 2014, reported in Table 1.
In order to obtain results comparable with those of Tigrisat mission analysis [18], the same values for the initial angular rates were set and reported below.
The control gain was set equal to so that the maximum available for the coils is not exceeded. The initial RAAN, orbit phase angle, and attitude are generated as random values within 0 and 360 deg.
The time behavior of the rotational kinetic energy is reported in Figure 1 from which it can be noticed that decreases by two orders of magnitude in less than 1 orbital period, within roughly 4500 seconds.
The following values for the maximum current () and torque () along the body axes were calculated as follows:
The algorithm is therefore suitable to produce the detumbling within the constraints defined for Tigrisat. An analysis of its performance is reported in Section 5.
4. Angular Rate Determination
As it is well known, the B-dot control is Lyapunov stable [14]. In virtue of it, under its effect, the magnitude of constantly decreases, approaching a value about zero. For a variety of missions, maintaining a residual spin of rotation could be desirable [19]. The capability of measuring the angular rates allows the detumbling control to be stopped when a given value of is matched. Measures of are here obtained using three consecutive samplings of from a three-axis magnetometer, fixed to the satellite rigid body and aligned with .
In accord to equations (4) and (6), the three vectors , , and are always mutually orthogonal and their directions, namely the unit length vectors, form a tri-orthogonal auxiliary reference frame, shown in Figure 2.
Since the values of are measured by a magnetometer fixed to the spacecraft, its direction in will change in time along with both the orbital position and the attitude of the rotating spacecraft. Considering now a short interval of time, corresponding to the inverse of the magnetometer sampling rate (), the variation of due to the spacecraft orbital motion can be neglected and, therefore, would rotate with , with its direction depending only on the spacecraft attitude motion. As detailed in Section 3, this assumption affects the measures of by a factor equal to the rate of rotation of the geomagnetic field vector, whose maximum value is equal to . Nevertheless, the proposed algorithm is targeted at spacecraft operating in LEO where 2n < 0.13 deg/sec. Therefore, the error introduced by this approximation is lower than the minimum measurement error achievable with the algorithm.
Using the superscript to indicate the generic sample of a fixed step time discretization, it will be proved that the following equation gives an accurate measurement of the spacecraft angular rate vector for the spacecraft controlled by the B-dot algorithm. where is the sampling rate and the time derivatives can be calculated through the upwind scheme
Equation (10) correlates the unknown angular rates of the spacecraft to the measurable rates of rotation of . Its validity is here proved for a spacecraft with a spherical mass distribution and principal moment of inertia . The results are then extended to the general case.
For a spherical spacecraft, (1) can be rearranged as follows:
Equation (12) shows that the magnetic control torque can only produce variations of along and, as proved hereafter, these variations correspond to rotations of the vector around .
According to (3) and (9), the vector always lies on the [] plane. Recalling that for adequately high, we can assume that depends only on the attitude motion of the spacecraft, and any rotation of around must be caused by a rotation of the spacecraft around the same axis (see Figure 3). Hence, a variation of along corresponds to the rate of rotation defined by (10).
For a spacecraft with a generic mass distribution, (1) must be considered in its original form and the direction of cannot be defined a priori. Nevertheless, also in this framework, the rate of rotation of can lead to angular rate determination. In fact, with respect to the case examined before, will also change along and . The components along these two directions are due to the shape of the tensor of inertia and to the nonlinear terms in the rhs of (1). When calculating the angular rates by means of (10), the and components introduce a measurement error that has the shape of a high frequency noise. Figure 4 shows the measurement error on the rotational kinetic energy (), defined as follows:
The measurement accuracy can be increased by compensating the effects of the nonlinear terms and filtering the results using a second-order Bessel low-pass filter. The compensation is obtained by adding to (10) the following one:
The results from (10) and (14) are then processed through a second-order Bessel low-pass filter that is characterized by the following transfer function: where is the filter cutoff frequency in Hz, whose value is unknown a priori and can be determined through simulations.
For the numerical analyses discussed in Section 5, the values of were selected in accord to the following criterion. First, the normalized (nondimensional) cutoff frequency () of the filter is defined as follows:
Then, a reference value for the normalized cutoff frequency (), for the axis corresponding to the mean principal moment of inertia (), can be determined through numerical simulations, assuming is equal to the maximum angular rate that can be expected at the deployment (this value is eventually provided by the deployer datasheet or can be guessed during mission analysis, depending on the geometry of the satellite and on the characteristics of the deployer).
Finally, the other values of can be calculated scaling by a factor proportional to
Angular rate determination was tested, along with the detumbling control discussed in Section 3, by means of numerical simulations in MATLAB. Both the angular rates obtained from the numerical integration of the attitude dynamics equations (solid line) and from the algorithm (dashed line) are reported in Figure 5.
The difference between the two, here regarded as the measurement error, can be seen in Figure 6, showing that the steady state error is limited within ±0.2 deg/sec.
Angular rate measurements can be particularly useful for missions in which different control on the three axes is required; this is for instance the case of Y-Thompson spin [20, 21]. An application of such a control was verified as well assuming a referenced angular rate along the y-axis is equal to 2.5 deg/sec.
The time behavior of the angular rates determined from both the attitude dynamics equations (solid line) and the algorithm (dashed line) are reported in Figure 7.
It can be noticed that both and are damped to zero, while only decreases down to the desired value. The measurement errors are shown in Figure 8.
The case study allowed validating the algorithm. Its performance and applicability on multiple scenarios with uncertainties are discussed in the next section.
5. Simulation with Unknown Deployment Conditions and Uncertainties
The target is to design an algorithm which is effective even when the deployment conditions are unknown or uncertain. This was proved simulating the detumbling for a variety of scenarios, in which the initial RAAN, orbit phase angle, attitude vector, angular rates, inclination, and altitude are unknown. Also, a 10% error on was assumed. In fact, the moments of inertia of a CubeSat are sometimes known with little accuracy, given the complicate and time-expensive procedure necessary for determining them. The mentioned constraints are hereafter summarized as follows:
The value of the control gain matching all the mentioned requirements was defined from the simulations, resulting in . This value allows a satisfactory detumbling for all the scenarios considered and was set constant for the sake of simplicity, although an optimal value for the gain might be determined [22]. The parameters of the low-pass filter, defined as in Section 4, are summarized in Table 2.
A total of 100 test cases were simulated for a time equivalent to 3 orbital periods. The angular rates and the measurement errors for all the cases are shown in Figures 9 and 10, proving that the algorithm is suitable to produce both detumbling and the angular rate determination regardless of the uncertainties on the deployment conditions and on the inertial properties of the CubeSat.
Once the robustness of the algorithm with respect to the parameters considered was verified, the effect of each single one was investigated. In order to compare the results from different test cases, two evaluation parameters are introduced, the detumbling time (), here defined as the time to reduce by two orders of magnitude the initial and the settling time (), namely the time at which the measurement error enters and remains within the error band ± 0.2 deg/sec.
We first focus on the effect of the deployment conditions, assuming therefore an exact knowledge of the tensor of inertia and a fixed value for the sampling rate, equal to 10 Hz. Simulations for these conditions were performed, and relevant results are reported in Tables 3 and 4.
The analysis of these data suggests that is influenced by the initial value of and . In particular, the variations in are more marked for higher values of the corresponding principal moment of inertia.
It is possible to notice as well that as the inclination of the orbital plane increases, the detumbling time decreases while the settling times are almost constant.
The influence of on the accuracy of the angular rate determination algorithm is now examined. According to (10) and (11), increasing the sampling rate corresponds to a reduction in the accuracy of the time derivatives calculated and, consequently, on that of . An example of this effect is shown in Figure 9, representing the measurement errors for 3 different values of (10 Hz, 8 Hz, and 1 Hz).
As shown in Figure 11, increasing the sampling rate from 1 to 10 Hz, the measurement errors grow and the effect seems to be more marked on the x-axis, the one with the minimum moment of inertia. Despite the lower accuracy, also for the highest value of , the measurement errors rapidly converge to the ±0.2 deg/sec error band.
Finally, the uncertainty on is taken into account. The measurement errors, calculated assuming exact knowledge of , a 10% error on and a 10% error on both and are represented in Figure 12. These results show that under- or overestimating the values of by one-tenth determines an increase of the mean measurement error in the early filtering phases, though not significantly affecting the performance of the system. This behavior can be better visualized in Figure 13, reporting the measurement errors for the initial 600 seconds of the detumbling. It is worth to outline that, within the range examined, an error on one moment of inertia only affects the measurement of the angular rate on the corresponding axis.
As stated at the beginning of the manuscript, the purpose of the present algorithm is allowing the detumbling and angular rate determination of a CubeSat based on the only data from a three-axis magnetometer. The results from the analyses presented in this section indicate that the algorithm is suitable for both producing the detumbling under unknown or uncertain conditions and measuring the angular rates in real time, with a measurement error rapidly converging to ±0.2 deg/sec. According to it, the algorithm can represent for CubeSats an adequate alternative to rate sensors and a suitable backup solution in case of their failure, eventually allowing the recovery of the mission.
Once the performance of the algorithm is verified, a comparative analysis can be performed, to characterize it with respect to similar solutions available in scientific literature, highlighting the novelty it introduces and identifying its limitations. We recall here that gyroless algorithms for angular rate estimation can be classified into two main classes, estimator-based, more accurate but computationally expensive, and deterministic, more efficient and easy to be implemented on a CubeSat onboard computer [3]. The proposed work belongs to this second class.
Deterministic gyroless algorithms typically require attitude knowledge, either directly measured or calculated, to be differentiated and used along with the attitude kinematic equations to compute the angular rates. Because of the noise introduced by the differentiation process, further data processing is necessary to improve the accuracy of the results and it can be obtained by means of either an estimator (usually a Kalman filter) or a passive one, such as the low-pass filter discussed in this work, significantly reducing the computational cost.
The use of a low-pass filter was already proposed by Bar-Itzhack, for processing the coarse angular rates calculated differentiating the attitude data determined from the measurements of the three-axis magnetometer and the sun sensor on the Rossi X-ray Timing Explorer satellite [23]. It can be noticed that the behavior of the estimated angular rates reported in the referenced paper, before and after the filtering, resembles those shown in the previous sections. Further development by Bar-Itzhack et al. [24], reporting rate estimation error graphs, allows a more detailed analysis. The plots related to the simple passive filter show that rate estimation errors are comparable to the measurement errors presented by the authors in this manuscript. Despite the similarity in terms of accuracy, it is worth noting that our results are obtained based on the only inputs of a single three-axis magnetometer and without measuring or computing the attitude.
The lack of dependence on attitude knowledge and the use of the only three-axis magnetometer are therefore the most peculiar aspects of the work presented in this manuscript. A similar magnetometer-only deterministic algorithm was originally proposed by Oshman and Dellus [25], for a LEO satellite actuated by momentum wheels, and later improved by Psiaki and Oshman [26]. The angular rates are here determined calculating the orthogonal component by means of an equation equivalent to (4) and the parallel component from the global solution of a least-squares problem. The same result in our algorithm is obtained by means of only elementary algebraic operations. This follows from the different attitude dynamics of a spacecraft controlled by magnetorquers, our case, instead of reaction wheels, as in [25, 26]. The actuation by means of a magnetic torque implies the geometrical constraints, discussed in Section 4, between the vectors , , and , finally resulting into (10).
6. Conclusions
The design and numerical validation of a magnetometer-only detumbling and angular rate determination algorithm was presented. A variety of deployment scenarios was numerically simulated and analyzed, showing that the proposed algorithm is suitable for yielding three-axis magnetic stabilization of a CubeSat deployed in LEO with unknown RAAN, orbit phase angle, attitude, inclination, and angular rate. Assuming a 3 U CubeSat as reference and fixing the maximum commanded magnetic dipole moment to 0.3 Am2, it was proved that the detumbling can be achieved within 2 orbital periods.
Angular rates can be determined in real time, based on three-axis magnetometer measurements. Data collected are processed, compensating the effects of nonlinearities, and filtered through a second-order low-pass filter, obtaining a final measurement error within the error band ± 0.2 deg/sec. The robustness of the algorithm, with respect to unknown deployment conditions, sampling rates ranging from 1 to 10 Hz, and uncertainties up to the 10% on the moments of inertia of the CubeSat are verified. Analyses showed that higher measurement errors correspond to higher sampling rates, producing an extension of the settling time. With respect to the moments of inertia, uncertainties on their value up to the 10% produce a slight increase of the measurement error in the early filtering phase but do not have relevant impact on the settling time.
Allowing angular rate determination, the algorithm can be conveniently used to monitor the detumbling phase, especially in those cases when a residual spin of rotation around an axis is desired. Therefore, by the algorithm, a convenient condition to switch the attitude control of the CubeSat from detumbling to fine pointing or attitude motion can be determined, eventually producing a considerable reduction in electric power consumption.
Nomenclature
| Bb: | Geomagnetic field in body-fixed coordinates |
| B i: | Geomagnetic field in GCI coordinates |
| Fb: | Body-fixed reference frame |
| Fi: | Inertial reference frame |
| fk: | Sampling rate |
| h: | Spacecraft altitude |
| [I]: | Identity matrix |
| I: | Electric current through the coil |
| i: | Inclination of the orbital plane |
| [J]: | Spacecraft tensor of inertia |
| Jr: | Reference principal moment of inertia |
| K: | Control gain |
| m: | Magnetic dipole moment |
| n: | Spacecraft orbital rate |
| T: | Orbital period |
| Tr: | Rotational kinetic energy |
| td: | Detumbling time |
| θ: | Spacecraft pitch angle |
| λ: | Cutoff frequency |
| τ: | Control torque |
| φ: | Spacecraft roll angle |
| ϕ0: | Orbit phase angle |
| ψ: | Spacecraft yaw angle |
| Ω: | RAAN |
| ω: | Spacecraft angular rate vector |
| : | Projection of ω orthogonal to Bb. |
Conflicts of Interest
The authors declare that they have no conflicts of interest.