Investigation of Instabilities in Detumbling Algorithms
IInvestigation of Instabilities in Detumbling Algorithms
Jeet YadavBirla Institute of Technology andScience, PilaniRajasthan, India - [email protected] Tushar GoyalBirla Institute of Technology andScience, PilaniRajasthan, India - [email protected]
Abstract—
Detumbling refers to the act of dampening the an-gular velocity of the satellite. This operation is of paramountimportance since it is virtually impossible to nominally performany other operation without some degree of attitude control.Common methods used to detumble satellites usually involvemagnetic actuation, paired with different types of sensors whichare used to provide angular velocity feedback.This paper presents the adverse effects of time-discretizationon the stability of two detumbling algorithms. An extensiveliterature review revealed that both algorithms achieve absolutestability for systems involving continuous feedback and output.However, the physical components involved impose limitationson the maximum frequency of the algorithm, thereby makingany such system inconceivable. This asserts the need to performa discrete-time stability analysis, as it is better suited to reflecton the actual implementation and dynamics of these algorithms.The paper starts with the current theory and views on thestability of these algorithms. The next sections describe thecontinuous and discrete-time stability analysis performed by theteam and the conclusions derived from it. Theoretical inves-tigation led to the discovery of multiple conditions on angularvelocity and operating frequencies of the hardware, for whichthe algorithms were unstable. These results were then verifiedthrough various simulations on MATLAB and Python3.6.7. Thepaper concludes with a discussion on the various instabilitiesposed by time-discretization and the conditions under which thedetumbling algorithm would be infeasible. T ABLE OF C ONTENTS
1. I
NTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2. D
ETUMBLING A LGORITHMS . . . . . . . . . . . . . . . . . . . . .
3. C
ONTINUOUS -T IME S TABILITY A NALYSIS . . . . . .
4. D
ISCRETE -T IME S TABILITY A NALYSIS . . . . . . . . .
5. S
IMULATIONS AND D ISCUSSIONS . . . . . . . . . . . . . . . .
6. C
ONCLUSIONS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A CKNOWLEDGMENTS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . R EFERENCES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B IOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1. I
NTRODUCTION
Detumbling refers to the act of reducing the angular velocityof the satellite to come under an acceptable predefined value.The state of high angular velocity can be induced afterdeployment from the launch vehicle, or naturally in orbit dueto disturbance torques. Detumbling is an essential task, asit would precede any other operation on the satellite whichneeds some degree of attitude control.The work presented in this paper was conducted in associ-ation with Team Anant, a group of undergraduate studentsworking to build a 3U CubeSat. Specific to the mission, lowangular velocity is desirable as an entry condition to various $31 . c (cid:13) pointing and tracking modes [1]: (1) Hyperspectral Imaging,(2) GS Tracking, (3) Sun Pointing.The team explored two different detumbling algorithms andcompared their efficiency, stability and power metric. Thesimulations were performed for a 3U Cubesat in a sun-synchronous, Low Earth Orbit (LEO) using MATLAB [2]and Python 3.6.7. Since detumbling is a high priority task,a thorough stability analysis was done which aimed to checkthe feasibility of detumbling the satellite from any giveninitial condition and chosen timestep. As reported in [3], theexpected rate of initial tumbling is 10 deg/s. However, im-perfect deployment or space debris might cause an unusuallyhigh tumbling rate, which the satellite should be prepared totackle appropriately. An example of this is the Swisscube[4], which had an initial tumbling rate of 200 deg/s andwas allowed to naturally detumble for a year, before normaloperations commenced.The paper aims at analysing the stability from two differentperspectives. While the continuous time analysis does hintat the behaviour of the controller, and clearly the reasonsfor its robustness, we also see that it fails to mimic the trueimplementation in some cases. The discrete time analysisshows how initial conditions and the parameters chosen forthe controller might affect the stability of the algorithms.The next section goes over the essential background infor-mation needed to perform the analysis. The theory and onboard implementation of two algorithms are discussed. Thissection is followed by a summary of the existing perspectiveson stability for the respective algorithms. The subsequentcontent shows the limitations of this perspective, and theadvantages in shifting to a discrete-time analysis, which accu-rately reflects the actual implementation of these algorithms.Instability criteria are derived and explained for simple cases.Simulations are also performed which demonstrate the valid-ity of the theoretical approaches used, and show the behaviorof these controllers for asymmetrical bodies. The paperconcludes with a discussion on the instability conditions andthe results from the various simulations performed.
2. D
ETUMBLING A LGORITHMS
This section describes the background of the selected de-tumbling algorithms and the dynamical equations used toanalyse the behaviour of the satellite. While the algorithmsdiffer in the choice of sensors, both use magnetic actuation todetumble the satellite. Thus, both algorithms aim to calculatean appropriate magnetic moment which must be produced soas to interact with the external magnetic field and reduce theangular velocity.Given the choice of magnetic actuation, the expressions forthe torque and the rotational dynamics equation are as follows1 a r X i v : . [ ee ss . S Y ] S e p ˙ ω BIB = ( I B ) − [ Γ B − ω BIB × ( I B ω BIB )] (1) Γ B = m B × b B (2)In the above equations, ω BIB is the angular velocity of thebody with respect to an inertial frame, represented in thebody frame. Γ B , m B and b B is the external torque (onlydue to magnetic interactions), magnetic moment, and externalmagnetic field, all represented in the body frame. For the sakeof convenience, we will drop the subscript and assume thatthe body frame representation is used except when indicatedotherwise.It must be noted that, while using magnetic actuation, we cannever have a component of torque along the instantaneouslocal magnetic field. This can be seen as a result of (2). How-ever, the lack of control of the angular velocity componentparallel to this field is taken care of due to the spatial andtemporal variation of the magnetic field in an orbit [6].The first control algorithm uses both angular velocity andmagnetic field information as feedback. The second algo-rithm (BDot) uses only the latter as input. The rate at whichthese values are received from the magnetometer and IMUwill determine the frequency at which we can run thesealgorithms. Discretized controller and considerations for itsimplementation are discussed later. Algorithm 1: ( ω × b ) The magnetic moment generated is perpendicular to the an-gular velocity and the local magnetic field vector. Angularvelocity can be split up into two components: a componentwhich is along the direction of the local magnetic field and acomponent which is normal to it. The magnetic moment forthe control law is calculated as follows [5] m = k c (cid:107) b (cid:107) ( ω × b ) (3)Here, k c is a scalar gain, ω is the angular velocity of thesatellite, and b is the local magnetic field. This particular se-lection of magnetic moment ensures that the torque producedis antiparallel to the angular velocity component normal tothe magnetic field.Feedback for this control law comes from the Inertial Mea-surement Unit (IMU), as well as the magnetometer. It is tobe noted that the use of onboard IMU will require calibrationand evaluation of bias and drift. Algorithm 2: ˙ b The Bdot control law calculates magnetic moment using therate of change of the magnetic field. It utilizes feedbackexclusively from the magnetometer. The control law takesthe following form [5]. m = − k c (cid:107) b (cid:107) ( ˙ b ) (4)In actual implementation, the rate of change of magnetic fieldis calculated by using a finite difference method as describedin (5) ˙ b = b k − b k − T s (5) Comparison
It can be shown that the control law for both algorithms aresimilar, and in fact identical under the following assumptions.1. ˙ b I ≈ ˙ b B is a continuous time derivative, and not calculatedusing (5).The magnitude of the effective rotation rate of the exter-nal magnetic field is of the order of the orbital rate ( − rad/s). This is several orders lower than the angular velocityencountered when the satellite is tumbling, and hence theassumption that it can be ignored is valid. This assumptionbreaks down when the angular velocity of the satellite is verysmall. In such a case this rotation rate is significant andhinders the controllability from that point forward. Therefore,while similar to the first algorithm, Bdot cannot completelydetumble the satellite.The second assumption states that Bdot should be calculatedin continuous time, and that the finite difference method usedto calculate the former is not always a reliable representationof the angular velocity. In practice, the magnetometers alonecannot give us Bdot in continuous time and is limited byits maximum frequency of operation. Therefore, a smalltimestep is desired so as to mimic a near-continuous calcu-lation. Further limitations of the finite difference method arediscussed in the following sections.Substituting the aforementioned assumptions in equation (6),we can show that that the magnetic moment induced by bothalgorithms are indeed equal (7). ˙ b B = A BI ˙ b I − ω BIB × b B (6) ˙ b B ≈ − ω × b B (7) Scalar Gain
The scalar gain k c (3) (4), to be used in the detumblingalgorithms on board, was selected after a literature reviewand various simulations. A possible candidate for the sameis shown below.The gain expression, proposed in [7] is based on analysing theclosed loop dynamics of the component of angular velocityperpendicular to the earths magnetic field. k c = 4 πT orb (1 + sinξ ) J min (8)Here, T orb represents the orbital time of the satellite, ξ represents the inclination of the satellite with respect to thegeomagnetic equatorial plane, and J min is the minimumprinciple moment of inertia for the satellite.A close look at the gain (8) and the constants present in thecontrol laws reflect independence of the angular accelerationfrom the magnitude of the magnetic field and relative size ofthe body. This results in effective control which is strictly afunction of the current angular velocity.
3. C
ONTINUOUS -T IME S TABILITY A NALYSIS
In this section, we approach the problem of stability analysisby assuming continuous time feedback and control. The2otational dynamics concerning angular velocity is as givenby (1). For torque determined by the control laws, theanalysis of the dynamical system is expected to show that theangular velocity will asymptotically converge to zero.In some cases, the behaviour of a system can be inferred byusing Lyapunovs second method [8]. This method relies onthe observation that asymptotic stability is very well linked tothe existence of a Lyapunovs function V ( x ) , which is definedby the following conditions1. V ( x ) > for all x (cid:54) = 0 ; V (0) = 0 d ( V ( x )) dt ≤ for all x The stability thus described implies that the equilibrium point x ∗ = 0 is stable and locally attractive. If condition two ischanged to a strict inequality, it is shown that the existenceof a Lyapunov function for a given system guarantees asymp-totic stability.From this perspective, we proceed to define a Lyapunovfunction candidate for our system: V ( ω ) = 12 ω T ( I B ) ω (9)where I B is a symmetric positive definite matrix describingthe moment of inertia of the body. In a physical sense, V ( ω ) is the rotational kinetic energy of the rigid body. Thestrictly positive definite nature of this function comes from itsquadratic form and the properties of I B .To apply the second theorem of Lyapunov, we take a timederivative of this function: ˙ V ( ω ) = 12 ( ˙ ω T I ω + ω T I ˙ ω )= ω T Γ (10)We can get an expression for the torque by subsituting (3) and(4) in (2). Furthermore,(7) shows that the control laws forboth algorithms can be assumed to produce the same torque,which is the following: Γ = − k c ( − ˆ b ˆ b T ) ω (11)Subsituting (11) in the expression for ˙ V , we get ˙ V = − k c ω T ( − ˆ b ˆ b T ) ω = − ω T Φ ω (12)Where Φ is a positive definite matrix. Therefore ˙ V isexpressed in a negative semidefinite quadratic form, and iszero only when w tends to zero or ω is parallel to b . Its phys-ical interpretation indicates that the rotational kinetic energyalways decreases, except for when the angular velocity is inthe direction of the instantaneous magnetic field, or is zero.This result is consistent with the fact that we can never havea component of torque along the magnetic field.Once the angular velocity is be reduced only to its componentalong the field, the control is effectively stopped. However,as mentioned before, the magnetic field is time variant.The change in the field would misalign it from the angularvelocity, and hence control will be resumed. Therefore, inall practical situations, asymptotic convergence of angularvelocity to zero is achieved due to the global variation of themagnetic field.
4. D
ISCRETE -T IME S TABILITY A NALYSIS
The practical implementation of these control laws is indiscrete time, and the frequency of operation is dependent onthe hardware employed. The rate at which these algorithmsrun will be dependent on the maximal frequency of the sensoras well as the characteristics of the processor and actuators.In using these algorithms in discrete time (assuming dutycycle unity), the magnetic moment calculated is constant inthe body frame, for a given time step ∆ t . This implies thatthe change in torque over ∆ t is due to the relative motionof the magnetic field with respect to the constant magneticmoment.The discretized dynamical equation (13) describes a non-linear non-autonomous system, where ˙ ω ( t ) describes theangular acceleration acting on the body, as evaluated usingEq. 1. ω k +1 = ω k + (cid:90) t +∆ tt ˙ ω ( t ) dt (13)In this section, we will explore three types of instabilities.The analysis follows reasonable assumptions in order to sim-plify the equations and build an intuitive sense for the sourceof the instability. The exact dynamics and results achieved bynumerical methods are shown in the next section. Type I
This instability arises when the control torque applied overthe timestep causes certain components of the angular ve-locity to flip and increase,hence increasing the net rotationalkinetic energy.For the sake of analysis, we isolate this type of instabilityfrom the others mentioned by assuming that the angularvelocity of the satellite is low. Therefore the torque in thetime period ∆ t can be taken as constant. We also assume thatthe control law being used is (3), so that the inaccuracy of thefinite difference method used in Bdot is not included.We expect that the aforementioned instability will arise wheneither the timestep ∆ t or the magnitude of torque is toolarge. Under certain conditions described below, the controllaws defined above are not asymptotically stable. To useLyapunovs second method on a discrete system, we need toreplace the time derivative with a difference [9]. So criterionfor instability now becomes ˙ V ( x ) ⇒ V ( x k +1 ) − V ( x k ) > (14)For the purposes of analysis, we take the body to be spheri-cally symmetric. The [ ω BIB × ( I B ω BIB )] term is zero, andwe have effectively decoupled the components of angularvelocity. Given that there is no control along the b axis, thedynamical behavior is restricted to the plane perpendicular tothis axis. The discretized equation is now simply: ω k +1 = ω k + Γ ∆ tI (15)where Γ is the torque and I is the moment of inertia of thespherical body.By using the Lyapunov candidate function (9) and the relation(15), we can evaluate (14). This will reflect on the systemdynamics of angular velocity with discrete-time control. V ( ω k +1 ) − V ( ω k ) = 12 ω Tk +1 ( I ) ω k +1 − ω Tk ( I ) ω k (16)3ubstituting (11) as Γ , we get the following condition forinstability k c I ∆ t > (17)The condition is derived by isolating it from other sources ofinstability. We can recognize (17) in another form: Γ I ∆ t > − ω ⊥ (18)This form supports our initial expectations. It is clearlyseen that the controller is unstable if the torque acting on acomponent of angular velocity, enables a difference of morethan twice the component itself over that period of time. Thecomponent effectively flips in direction and keeps growing.This would occur in every iteration (when m is refreshed), andangular velocity would increase. This behaviour will persisttill the change in torque over the timestep is too large to beignored.This condition demonstrates the dependence of stability onthe timestep and gain chosen. Type II
This type is specific to the Bdot controller described in (4).The source of the instability lies in the limitations of usingmagnetometer readings to get a sense of the angular velocity.In the discretized controller, ˙ b in (7) is replaced by the finitedifference method shown in (5). From the equations, wecan see that the difference vector b k +1 − b k should give usinformation about the component of angular velocity perpen-dicular to it. Henceforth we shall refer to this component as ω ⊥ .The limitations occur when b rotates more than π radians.The discretized version of (7) will interpret that the smallerangle between b k +1 and b k is a result of ω ⊥ . Hence, whenthe angle θ is more than π radians, the algorithm assumes π − θ is the angle b rotated by. Therefore, the ω ⊥ is falselyassumed to be in the opposite direction, and the controllerwill contribute a positive angular acceleration to the satellite.This acceleration will continue until the controller reaches itsequilibrium point ω ∗⊥ = 2 π rad/s.The aforementioned limitation can be empirically illustratedby looking at the dependence of the induced angular velocity(at the start of the interval) on angular velocity. We willapproach this by taking the case of a spherically symmetricbody and analysing its behaviour. First, we assume a rotationabout the z-axis only. It should be noted that since the bodyrotates with ω with respect to the inertial frame, b will rotatewith − ω with respect to the body. Given a measurement b k − , the subsequent measurement after ∆ t is: b k = (cid:34) cos ( ω z ∆ t ) sin ( ω z ∆ t ) 0 − sin ( ω z ∆ t ) cos ( ω z ∆ t ) 00 0 1 (cid:35) b k − (19)By using (4) and (5), it can be shown that the expression forangular acceleration about the z-axis is ˙ ω z ( t k ) = − C sin ( ω z ∆ t ) where C = k c ( b x + b y )∆ t I (cid:107) b (cid:107) (20) Figure 1 represents the phase portrait with stable and unstableequilibrium points at even and odd multiples of π respec-tively, as described by (20). This shows us that in orderfor the satellite to start detumbling, and effectively reach theequilibrium point ω ≈ , the following criterion will have tobe fulfilled: ∆ t < π (cid:107) ω (cid:107) (21)It is worth mentioning that (21) is basically an expression of . π ∆ t π ∆ t . π ∆ t π ∆ t . π ∆ t π ∆ t . π ∆ t π ∆ t − kk Angular Velocity ( ω ) ˙ ω z ( t k ) Figure 1 . Stability and Instability Pointsthe Nyquist theorem. This theorem determines the maximumsampling time ∆ t needed in order to represent a signalwithout aliasing [10]. Hence, this maximum is half the timeit takes for b to rotate once.In conclusion, during the process of implementing BDot,attention also needs to be given to the maximum allowedangular velocity that the system can effectively control. Type III
This type of instability applies to both the algorithms. Itsbasis lies in the fact that the torque provided by the actuatorschanges in the timestep ∆ t .Due to the discrete nature of the update, magnetic moment isfixed for a given ∆ t . The rotation of the external magneticfield with respect to the induced magnetic moment, as seenin the body frame, results in a time varying torque. For de-tumbling to proceed, the integral of torque over this timestepshould be negative, so that there is a net deacceleration.To build an understanding of this effect, consider the body tobe spherically symmetric. As mentioned before, this restrictsthe dynamics to a plane and helps simplify the analysis. Thesystem of differential equations for this system, for the initialconditions shown in Figure 2, are as follows: ˙ θ = ω ˙ ω = − k c ω i I cos ( θ ) (22) ω and ω i are the angular velocity and the inital angularvelocity. The latter is a constant, while the former changesover the timestep. I is the moment of Inertia of the sphericallysymmetric body.As per the body frame representation shown in Figure 2, θ isthe angle by which magnetic field has rotated from its initialvalue b init to b . The magnetic moment m is fixed in the bodyframe.To discover our instability criterion, we now need to get anexpression for ω ( t ) . If this is smaller than the inital angularvelocity, then the satellite will detumble. However, gettingan exact expression for ω ( t ) is out of the question due to thenon- integrability of the system of differential equations (22).4 ωb initial b Oθ Figure 2 . Change in Magnetic Field, in Body FrameWe now assume that the time varying torque has a negligibleeffect on θ , as compared to the effect due to a high angularvelocity. The benefit of incorporating this assumption istwofold. Firstly, it helps us isolate this type of instability fromthe Type I, since condition (17) will never be fulfilled. Sec-ondly, this assumption can be used to make (22) integrable.While it should be noted that this assumption might not bevalid for all cases, it will help us look at the behaviour of thesystem during this instability.With this, we get θ = ω i t . Noting that the initial angularvelocity is fixed, the equations are solved ˙ ω = − k c ω i I cos ( ω i t ) ω ( t ) = − k c sin ( ω i t ) + ω i (23)The above equation illustrates the dependence of ω ( t ) on bothtimestep and initial angular velocity. Replacing t = ∆ t , wenow proceed to express the criterion for stability: ω ( t ) − ω i = ∆ ω = − k c sin ( ω i ∆ t ) < (24)It can be shown that, to satisfy the inequality in (22), thefollowing criterion will have to be fulfilled: ∆ t < π (cid:107) ω (cid:107) (25)After the theoretical approach given above, (22) is numeri-cally integrated and it is seen that both results match. Thefollowing graphs show the simulations performed by theteam, in order to validate the derived instability conditions.Figure 3 and 4 show the dependence of the change in an-gular velocity after one iteration, on both the initial angularvelocity, and the time step. The graphs seems to be ofthe same nature, which meets the expectations built using(23). Thus we can see that the governing factor for theinstability is actually the product ω i × ∆ t , which representsthe approximate angle by which the magnetic field has rotatedwith respect to the magnetorquers. This can be interpretedby paying attention to the fact that the torque producedstarts supporting the magnetic field after b enters the secondquadrant. This would then explain the points of zero change,where the same acceleration and decceleration cancel eachother out after rotations of nπ radians. Figure 3 . Change in angular velocity vs initial angularvelocity ( ∆ ω vs ω i ) Figure 4 . Change in angular velocity vs time step ( ∆ ω vs ∆ t )Keeping in mind the previously used assumptions, we nowdemonstrate that behaviour of the system is similar even withtorques of high magnitude. Figure 5 shows how the change inangular velocity is affected by increasing the coefficient of thealgorithm. We can see that both the curves follow a similartrend, with an apparent shift between them. Keeping in mindthat the torque can no longer be neglected to determine theangle, it is evident that a higher initial angular velocity isneeded to rotate by the same angle, for the case of the highercoefficient. It is to noted that the vertical axis have beenrelatively-normalized to explain the apparent shift.5 igure 5 . Change in angular velocity vs initial angularvelocity; for different Coefficients
5. S
IMULATIONS AND D ISCUSSIONS
This section consists of the various simulations performedto analyze and validate the various instability conditions, asmentioned in Section 4. It is to be noted that all the graphsshown here are for Algorithm 1 2, unless specified otherwise.We perfomed similar study for Algorithm 2 2 as well, and theresults were similar to the ones shown below.
Figure 6 . Angular velocity vs time step ( ω vs t )The figure 6 shows the variation of angular velocity with time,for different initial conditions. We can see that for variousangular velocities, the equilibrium point is different. Caseswith velocities between − π and π converge to zero, whilethose outside this range, converge to the nearest even multipleof π .The graph above is simulated for a fixed step time( ∆ t ) andcoefficient ( K c ). This indicates that any combination of theseparameters will have certain unstable points (not convergingto zero) , which have to be kept in mind while designing thealgorithms. Figure 7 . Average angular acceleration vs time step ( ˙ ω av vs ∆ t )Figure 7 shows the variation of the average acceleration ofthe body, for different step time used for implementation ofthe algorithm. It can be clearly seen that a lower step timeresults in better, faster desaturation of the angular velocity.This underscores the performance of the system in continuousdomain vs that in the discrete domain.The initial part of the trajectory follows our expectations andshows that the efficiency of the algorithm decreases, as wemove away from the continuous domain. After a point, themagnetorquer effect described in Subsection 4 takes over, andvarying amounts of acceleration and decceleration results inan oscillating curve, as described in (23) Figure 8 . Angular velocity vs time ( ω vs t )Figure 8 shows the variation of angular velocity over time,while using these algorithms for multiple time steps. Thisshows that the conclusions about the performance, drawnfrom Figure 7, holds true for longer periods of time as well.6 igure 9 . Average angular acceleration vs coefficient vsinitial angular velocity; Symmetric Body ( ˙ ω av vs K c vs ω i )Figure 9 is a result of the instability condition described inSubsection 4. We can see that for a symmetric body, with ∆ t/I = 1 , the controller becomes unstable at K c = 2 . Thevalley in the surface represents the point at which the angularvelocity is only along the magnetic field, and hence, at itsminimum.An important point to be noted is that while we assumed lowinitial angular velocity to derive the instability condition, it isclear that the conditions holds for arbitrarily high initial an-gular velocities as well. It is demonstrated that the instabilitycondition is independent of the angular velocity, as shown in(17). Figure 10 . Average angular acceleration vs step time vscoefficient; Anant Satellite ( ˙ ω av vs ∆ t vs K c )Simulating with a complex body, we numerically integrateover a long period of time. The variation of average angularacceleration with change in timestep and coefficient is shownin the surface of Figure 10. The results show little deviationfrom the derived and simulated case of a symmetrical body. Figure 11 . Average angular acceleration vs initial angularvelocity vs step time; Symmetric Body ( ˙ ω av vs ω i vs ∆ t )Figure 11 shows the change in the average angular accelera-tion, with varying step time and initial angular velocity. If wetry and look at the two components of the graph separately,we can make out two previously-established conclusions.For a particular high angular velocity, the angular accelera-tion changes with step time as shown in Figure 4. A similarcharacteristic can be seen here as well. For a specific steptime, the angular acceleration changes with the initial angularvelocity as shown in Figure 3. The same feature of the curvecan be seen in the surface plot as well.In fact, the ripples created in the surface re-trace the curve ω i × ∆ t = constant , which is exactly the result we estab-lished in Subsection 4.
6. C
ONCLUSIONS
The paper goes over multiple sources of instability in de-tumbling algorithms and derives criterion for the same. Theanalysis is done for the two algorithms introduced in thesecond section of this paper. The work was motivated by aneed to perform a discrete-time stability analysis and estab-lish checkpoints before entering detumbling mode.The first type of instability applies to both algorithms and isexplicitly dependent on the coefficient and time-step used.The criteria for the same can be evaluated on ground, pre-launch. The second type applies to only Bdot and relies onangular velocity and timestep chosen. This comes from thelimitations of the magnetometers and could be tackled bysampling the magnetic field more frequently. The third typeof instability applies to both algorithms and also depends onthe initial angular velocity and timestep chosen.After an initial comparative analysis, BDot was selected asthe detumbling algorithm for the satellite [11]. The instabilitycriterion for this algorithm was then determined. The firsttype of instability depends on the parameters of the controller,as decided in the design phase. It is seen that ∆ t = 1 ,along with the optimized gain proposed by [7] successfullyavoids the pitfalls of this type. Given the timestep mentionedabove, the second and third instability can only be inducedby variations in the initial conditions. Using the derivedcriteria (21), (25), and giving some margin of error due to theasymmetric body of the satellite, we decided on a condition7or entry into the detumbling mode. If angular velocity ishigher than a predefined value, it will not be useful for thesatellite to enter into this mode. It should be allowed tonaturally detumble through drag in such a scenario.The discussed instabilities underscore the need to ensure thatthe control algorithms designed by any satellite manufacturerare robust. The discussions and results of this paper will beof essential assistance in ensuring that robustness. A CKNOWLEDGMENTS
The authors would like to thank the following:1. Team Anant, the Student Satellite Team of BITS Pilani,for providing the motivation and resources for the work donein this paper.2. Dr. Kaushar Vaidya, Vishnu Katkoori, and Jivat Neet Kaurfor their constant support and technical expertise in perusingresearch for this paper.3. The administration of BITS Pilani and Indian Space Re-search Organization (ISRO) for giving us an opportunity towork on building a satellite. R EFERENCES [1] Rutwik N. Jain, et al.
Modes of Operations of a 3U Cube-Sat with a Hyperspectral Imaging Payload . InternationalAstronautical Congress 2018, Bremen, Germany[2] Tushar Goyal, et al.
Simulator for Functional Verifi-cation and Validation of a Nano-satellite . 40th IEEEAerospace Conference 2019, Big Sky, Montana, USA.DOI: 10.1109/AERO.2019.8741886[3] Young-Keun Chang, et al.
Rapid Initial DetumblingStrategy for Micro/Nanosatellites Using a Pitch Bias Mo-mentum System . Transactions of The Japan Society forAeronautical and Space Sciences, 50(167):63 69, May2007.[4] Ecole polytechnique federale de Lausanne.SwissCube Project - The first Swiss Satellite.http://swisscube.epfl.ch/, September 2014.[5] F. Markley and J. Crassidis,
Fundamentals of spacecraftattitude determination and control
Tumbling of arigid spacecraft . Journal of Guidance, Control, and Dy-namics, 35(4):13261334, 2012. doi: 10.2514/1.53074.[8]
Lyapunov Stability Theory ,http://web.tuat.ac.jp/ ∼ gvlab/ronbun/ReadingGroupControl/mls-lyap.pdf[9] Lyapunov Theory for Discrete Time Systems ,http://automatica.dei.unipd.it/tl files/utenti2/bof/Papers/NoteDiscreteLyapunov.pdf[10] Alan V. Oppenheim and Alan S. Willsky,
Signals andSystems 2nd ed , Prentice-Hall 1997[11] Vishnu P. Katkoori, et al.
Simulation and Selection ofDetumbling Algorithms for 3U CubeSat . InternationalAstronautical Congress 2019, Washington DC, USA B IOGRAPHY [ Jeet Yadav , a student of Birla Insti-tute of Technology and Science, Pilani,India is pursuing Dual Degree in M.Sc.Physics and B.E.(Hons.) in Electron-ics and Electrical Engineering. He hasbeen a part of Team Anant since May2018, contributing to the developmentof Attitude Determination and ControlSubsystem (ADCS). He has been sub-system lead of ADCS since November2018, managing and working on various projects related toADCS and system integration. His interests lie in TheoreticalPhysics and he wishes to pursue research in the same.
Tushar Goyal received his B.E.(Hons.)in Electrical and Electronics Engineer-ing, with an unofficial Minor in Physics,from Birla Institute of Technology andScience, Pilani, India in 2019. He wasactive member of Team Anant, from Jan-uary 2016 to May 2019. He workedon the development of the Attitude Con-trol System overall system integration ofthe satellite. He is currently a SpaceDynamics and Systems Engineer at Astrome Technologies,India, working in various domains related to Space Tech.His interests lie in Astrodynamics & Mission Analysis andhe aspires to pursue doctoral research in the future.received his B.E.(Hons.)in Electrical and Electronics Engineer-ing, with an unofficial Minor in Physics,from Birla Institute of Technology andScience, Pilani, India in 2019. He wasactive member of Team Anant, from Jan-uary 2016 to May 2019. He workedon the development of the Attitude Con-trol System overall system integration ofthe satellite. He is currently a SpaceDynamics and Systems Engineer at Astrome Technologies,India, working in various domains related to Space Tech.His interests lie in Astrodynamics & Mission Analysis andhe aspires to pursue doctoral research in the future.