A Variable Sampling-Time Method for Elliptical Orbit Motion Prediction in Nanosatellites

Nanosatellite missions may contain payloads for high pointing accuracy such as laser communication systems for crosslinks or astronomical observations. Therefore, the satellite requires a precise orbital position and orientation determination in order to point the scientific instrument to the desired target. In this work, an elliptical rotation method based on quaternion representation is presented. The proposed method allows determining the future position of a satellite around its orbit. Furthermore, in Low Earth Orbits (LEO) with an eccentricity larger than zero, the distance between the satellite and the Earth is changing over the time, increasing the satellite velocity in the perigee region compared to the apogee, due to the gravity forces. The elliptical rotation method and the orbital current position are deduced, considering a variable sampling-time as a function of the eccentricity and the orbital current position. The proposed algorithm can increase the position accuracy four times compared to fixed sampling time along the satellite orbit.


I. INTRODUCTION
In an artificial satellite, the Attitude and Orbit Determination and Control Subsystem (AOCS) computes the current satellite orientation and its current position vector around its orbit. There are satellite payloads that may require a pointing accuracy of millidegrees such as laser communication in links or astronomical observations. In Low Earth Orbit (LEO), the Earth's magnetic field and the magnetic dipole moment of a nanosatellite produce a disturbance torque in the order of 10 −4 to 10 −5 Nm. For a nanosatellite, the magnetic disturbance [1] cannot be neglected to design a pointing control system in the order of millidegrees, since the nanosatellite mass is of 1 to 5 kg [6]. Therefore, disturbances can introduce a pointing error from 1 to 10 millidegrees.
In [2], a Fuzzy Variable Structure Control (FVSC) is proposed. The FVSC achieves a stable response considering acquisition delays due to sensors response time. The FVSC was designed for small satellites with an accuracy of 2 • . However, if higher accuracies are required, it would be necessary to increase the set of rules and this will affect the computational complexity. On the other hand, The associate editor coordinating the review of this manuscript and approving it for publication was Venkata Ratnam Devanaboyina .
PD and PID controllers are frequently cited in the literature [3]- [5] since they have high stability and fast response time. The PD and PID controllers alone can reach an accuracy up to 10 −2 degrees. In [7], an Extended Kalman Filter (EKF) is proposed, the authors conclude that the EKF requires more tunning to obtain high accuracy for attitude determination. A nonlinear stochastic filter is presented in [9], [10], the authors propose an attitude estimator based on Special Orthogonal Group SO (3). This algorithm computes the bias and noise produced by the sensors in order to minimize the attitude error. However, the complexity of this algorithm can be increased and rather large for a nanosatellite.
The principal motivation of this work is to propose an elliptical rotation model to estimate the future position of the nanosatellite. The proposed algorithm can be implemented aboard a nanosatellite in order to determine the future disturbance vectors which affect the attitude accuracy.
On the other hand, the rotation representation, formulated with Euler Angles, is susceptible to gimbal lock, since two of the three degrees of freedom can be along the same axis [12]. Since the quaternion representation does not present singularities, as in the case of Euler angles, the proposed model is based on quaternions in order to estimate the nanosatellite trajectory at any orbital inclination angle. With this goal, VOLUME 9, 2021 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ a model based on variable sampling-time is presented. The proposed model can adjust the response time of the AOCS system, according to the angular orbital motion, which adapts the sampling density where the satellite velocity is higher during its orbit.

II. REFERENCE FRAME FOR SATELLITES
A quaternion q ∈ R 4 consists of four elements [12]; the first three include the axis rotation vector (u) and the fourth element represents the rotation angle (ϕ). The quaternion is described as: where u ∈ R 3 . The rotation matrix A(q) ∈ R 3×3 and it is expressed as (2), shown at the bottom of the page. In this way, the rotation matrix A(q) can describe the satellite motion around its orbit as shown in Figure 1. On the other hand, the Keplerian orbit elements are the main parameters to mathematically describe the orbit propagation [8]. These elements are described in a two-line set, known as Two Line Elements (TLE), which can be updated when the satellite passes over an Earth station. The Keplerian elements are scalars and they are described as:  In Figure 2, the geometric elements of an orbit are represented on the Earth-Centred Inertial Frame (ECI). The point in the satellite orbit which is closest to the Earth is called the perigee and the apogee is the point farthest from the Earth.
The ascending node is the point of the orbit where the satellite crosses the equatorial plane from South to North and the descending node, where the satellite moves from North to South. The line of nodes is a straight line that crosses the ascending and the descending nodes. That is, it is the line that intersects the equatorial and the orbital planes. The argument of perigee is the angle between the ascending node vector and the perigee vector. The major semi-axis is the major distance between the satellite and the Earth center. The eccentricity describes the deviation of an orbit from circularity and the inclination is the angle between the orbital plane and the equatorial plane [8].
To determine the position of the satellite around the orbit, the mean anomaly (M ) is computed as: where E ∈ R denotes the eccentric anomaly. Moreover, it is required to compute, as a function of time (referred to as propagation), the longitude of the ascending node, the argument of perigee and the mean anomaly [13], respectively, given as: where η = µ a 3 is the mean angular velocity R e = 6, 378, 139 meters, is the Earth's radius. µ = Gm e = 3.986 × 10 14 m 3 s 2 , G is the universal gravitational constant m e is the Earth's mass. The zonal harmonics of the Earth (J n ) are computed taking into account the motion of the satellites around the Earth. The parameter J 2 is the strongest perturbation due to the Earth's shape: where J n ∈ R for all n ∈ N. In orbital propagators, the eccentric anomaly (E) and the true anomaly (M T ) are necessary to determine the position and velocity vector of the satellite. The eccentric anomaly (E) is given as: where the initial condition is E 0 = M . On the other hand, the true anomaly, M T (t) ∈ R, can be expressed as: The distance between the center of Earth and the satellite is defined as: where r c (t) ∈ R. In an orbit propagation based with fixed sampling time [13], the position vector and the velocity vector in the orbital frame are, Finally, the vectors o (t) andȯ (t) can be transformed to the ECI frame using the rotation matrices [11].
where r (t) ,ṙ(t) ∈ R 3 and they describe the position vector and the velocity vector of the satellite, respectively.

III. ROTATION ALGORITHM BASED ON VARIABLE SAMPLING-TIME
In elliptical orbits, the satellite varies its distance with respect to the Earth over time. Since the perigee is the closest point to the Earth, the satellite speed increases. In the AOCS subsystem, based on a fixed sampling-time, the angle difference is larger in the perigee in comparison to the apogee region. This angle difference produces a delay in the estimation of the future position, as shown in Figure 3. To correct the delay in the position estimation when the satellite is around the perigee, a predictive rotation model based on variable sampling-time (t s ), described in Equation (14), is proposed. In this way, the future position is estimated according to its proximity to the perigee and the eccentricity of its orbit, through the following expression of the sampling-time: where e is the eccentricity of the orbit, θ ∈ R is the angular position around the satellite orbit. If θ is equal to 2π, then the satellite is in the apogee region whereas the satellite is VOLUME 9, 2021 in the perigee region when θ is equal to zero and finally, t r describes the response time of the hardware. In order to calculate the parameter t r , we propose a normalized response time (t n ) as a function of the normalized sampling time of the sensors (t s ) and actuators (t a ), where K , L ∈ N are the total numbers of sensors and actuators, respectively. Thus, the parameter t ms is the normalized response time of sensors and t ma is the normalized response time of actuators. In fuzzy logic, the intersection can be related to the minimum function whereas the union can be associated with the maximum function [14]. In this way, from Equation (15) to (17) can be substituted as follows: and thus, t r is the denormalized of the response time (t n ). Later, the ascending node vector [x a , y a ] and the perigee vector x p , y p are required in order to determine the angular position ( θ ). These vectors are formulated according to the next equations: x a = cos ( ) , where x p, y p , x a , y a ∈ R and they are represented on the ECI frame. The vector normal to the equatorial plane can be expressed in quaternion representation as follows: where q o ∈ R 4 defines the quaternion of the orbit according to the inclination parameter (i). Keeping in mind Equation (1), the ascending node vector [x a , y a ] substitutes the axis rotation vector and the inclination (i) expresses the angle rotation. Through the quaternion, the perigee vector expressed on the tree axes is defined by: The vector r p can be updated according to the TLE data (by receiving the data set from the Earth Station). In this way, the current state vector of the satellite (r n ), obtained by the GPS sensor, and the perigee vector (r p ) can determine the angular position ( θ ) of the satellite around its orbit: where r p ,r n ∈ R 3 and n is the current sample.
In order to demonstrate that k = 1+e 1−e 2 −1, Equation (14) is formulated as: Assumption 1: Let k ≥ 0 ∈ R, which denotes the smooth response according to the angular velocity of the satellite. The parameter k allows adjusting the sampling time of the satellite around its orbit.
If an eccentricity is not equal to zero and the angular position ( θ ) is equal to π , then the satellite is in the apogee and the satellite reaches the minimum angular speed (ω min ). Substituting θ = π and t s = 2π/ω min into Equation (28), yields the following expression: On the other hand, the shortest possible time t r is equal to 2π /ω max , when the satellite achieves the maximum angular speed (ω max ). Therefore, the relation between the eccentricity, the minimum and maximum angular speed is formulated as: If the eccentricity (e) is zero, then ω max = ω min . On the other hand, if e is not zero, then ω max > ω min . Therefore, ω max ω min ≥1 and Equation (30) satisfies the condition that the term k has to be equal to or larger than zero.
Assumption 2: The parameter k can be formulated according to the angular momentum law, which describes that the angular momentum on the perigee region (r min v max ) is equal to the product of distance and velocity at the apogee, where r min , r max ∈ R denote the minimum and maximum distance between the satellite and the Earth, respectively. The maximum velocity is expressed by v max and the minimum velocity is represented by v min . The tangential velocities can be related to the angular velocities as follows: v max = r min ω max , If the orbit is circular, the eccentricity is equal to zero. In this case, the angular speed in the apogee has to be equal in the perigee. Hence, when e = 0, the Equation (36) is Finally, when the eccentricity is larger than zero, the ω max /ω min ratio is larger than one. In this way, Equation (36) can be substituted into Equation (30), giving the parameter k, as follows: where Equation (38) is related to the expression of the variable sampling-time, described in (28). Based on this proposed sampling-time, we can minimize the angle difference of the nanosatellite trajectory due to its velocity change in elliptical orbits. The rotation axis (u) of the orbit plane is computed to predict the path orbit, with the next expression: Later, the angular variation ( ) is computed between the current position vector (r n ), obtained by a GPS sensor, and the last position vector (r n−1 ) of the satellite in the Earth-Centred-Earth-fixed frame (ECEF), where ∈ R and r n ∈ R 3 for all n ∈ N and u z is the z component of the vector r n . When the orbit is not polar, the angular variation ( ) can be determined using the first two conditions of Equation (40). On the other hand, if the orbit is polar, the angular variation ( ) is computed through the last four conditions. Substituting Equations (39) and (40) into (1), the quaternion of the satellite, which describes the rotation of the position vector around the orbital plane according to the angle rate ( ), it is obtained as: where q s ∈ R 4 denotes the quaternion of the satellite rotation around its orbit. In the literature, the unit-quaternion is described as non-unique [11]. The rotation matrix can be obtained by two different unit-quaternion vectors, that is, A (q) = A(−q). To analyze the non-unique property, the Equation (41) multiplied by -1 is formulated as: The quaternion −q s describes the satellite rotation, around its orbit, considering opposite angle rotation with opposite axis rotation. For this reason, Equation (40) determines the rotation angle according to the direction of the axis rotation.
Based on the previous process, the future position vector is obtained by the following equation: wherer n is the position vector described by the GPS sensor and A (q s ) is the rotation matrix in quaternion form. When the satellite is in orbit, the GPS sensor can be affected by a bias (β). The bias can be produced by the measurement of the GPS and unknown parameters (due to environmental disturbances). In this way, Equation (43) is formulated as: where β ∈R 3 . The angular error ( ) between the current position (r n ) of the GPS sensor and the next estimated position vector (r n+1 ) can be determined for each axis, as follows: where = x , y , z , and ∈ R 3 . Finally, β can be computed as the mean of the angular variation error ( ), whose components are given as: where N ∈R is the number of total samples in one orbit. The components β x , β y and β z can be experimentally determined during one orbit. In this way, the bias can be subtracted to increase the angular accuracy.

IV. RESULTS
To test our proposal, the angular accuracy and the number of total cycles are analyzed for one satellite orbit. In simulations shown from Figure 4 to Figure 6, an eccentricity of 0.001 is considered, which is a typical eccentricity for CubeSats whereas the inclination is 0 • . The semi-major axis is 6,671 km. Figure 4 shows the angular variation, during one orbit, between the proposed elliptical rotation with variable sampling-time and an orbit propagation based on fixed sampling-time, described in Section II. The response of the variable sampling-time has less standard deviation in comparison to the fixed sampling-time from 1s to 0.25s, increasing the accuracy of the angular orbital motion.
In Figure 4, the angular accuracy of the variable sampling-time and the fixed sampling-time are compared considering different initial sampling-time t r from 0.1 seconds to 1 second. In other words, t r is the sampling-time required to obtain a valid response of the hardware, which is determined by the time response of the sensors and actuators, expressed in Equation (20). In Figure 4, the mean of each angular variation is subtracted, respectively, in order to compare the standard deviation of the different fixed sampling-time and the variable sampling-time. Figure 5 shows that the variable sampling-time maintains constant accuracy of the angular variation, considering different performance times for hardware; whereas the fixed sampling time causes that the accuracy of the angular variation varies linearly with respect to the performance time for hardware.
In order to estimate the total number of cycles required for an attitude control, the cyles for the proposed algorithm, during one orbit, are multiplied by the complexity of the possible subprocesses according to the AOCS system. In this way, the number of total cycles is expressed as: where C total is the number of total cycles of the AOCS system, A is the complexity of subprocesses of an attitude control algorithm and C t is the total cycles of the elliptical method based on variable sampling-time or fixed sampling-time. Figure 6 shows the total cycles for one orbit. The simulation is done to evaluate the required processing cycles of the variable sampling-time and the fixed sampling-time considering one orbit. We observe that the cycles based on a variable sampling-time do not depend on the sampling-time for the hardware. On the other hand, the fixed sampling-time shows a relation as a negative exponential function between the cycles and sampling-time t r . In Figure 7, the accuracy based on variable sampling-time tends to maintain the mean of the accuracy for altitudes from 300km to 600km because of the parameter k adjusts the sampling-time based on the angular momentum law and the eccentricity. The order of the accuracy varies according to the eccentricities from 0.0005 to 0.01; whereas the accuracy based on fixed sampling-time decreases whether the satellite is near to the Earth due to the velocity rate of the satellite for different altitudes.

V. CONCLUSION
In this paper, an elliptical method for satellite position estimation is proposed. The algorithm with variable sampling-time increases the angular orbital accuracy up to 4.4 times more than a traditional propagator with a fixed sampling-time. For the simulations, typical CubeSats orbit eccentricities were considered, from 0.0005 to 0.01, and altitudes from 300km to 600km. Moreover, when the satellite is in orbit, the eccentricity parameter can vary due to the atmospheric drag, and thus, the satellite changes its velocity in the perigee than in the perigee. The elliptical rotation method can adjust the variable sampling-time according to the eccentricity. Therefore, the proposed method can be implemented for nanosatellite AOCS systems with high accuracy requirements.
The proposed method also adjusts the time response of the AOCS system based on the angular orbital motion as a function of the eccentricity and according to the TLE over time. In this way, the elliptical rotation method based on a variable sampling-time keeps the factor of the total cycles (C total ) independent of the sampling time for hardware (t r ). Furthermore, the algorithm based on variable sampling-time can be implemented on low power consumption microprocessors, which implies a low cost and small size for the nanosatellite design.
A future work is to implement the proposed method based on variable sampling-time in hardware in order to evaluate the experimental number of cycles and to be evaluated in real environmental disturbances.