Simulation of the Attitude Control System on a CubeSat U1 Using Magneto-Torquers

Gustavo Mendoza-Torres, Yolanda Elinor Bravo-Garcia, Edmundo Reynoso-Lara, Hilario Pajaro-Adrian

Facultad de Electrónica, Benemérita Universidad Autónoma de Puebla, Puebla, México.

DOI: 10.4236/aast.2025.101001                

Abstract

In this work, we simulate the proportional derivative (PD) control of the control system for the Attitude of the CubeSat U1 based on a set of 3 orthogonal magnetic torsion coils called Magneto-Torquer and analyze its functionality. The simulation of the PD control of the Attitude system is based on the Propat simulator programmed in MatLab, created by Dr. Valdemir Carrara. We obtain the data that will allow us to position the X face of the CubeSat U1, pointing towards Puebla city. The electromagnetic theory that allows us to generate the magnetic torque in the Magneto-Torquer is described as the variables that will be controlled in its application to control the CubeSat. We consider the radius of the Earth to be 6,371,000 meters. In particular, it is proposed to simulate with an orbit altitude equal to 300,000 meters with an eccentricity of 0.0002, with polar orbit, using a Magneto-Torquer to generate the movements necessary to locate the CubeSat, based on the simulation. For the simulation, the real conditions of outer space and the Earth’s magnetic field at the location of Puebla are considered, and the values for positioning the X face of the CubeSat towards this region are obtained. Confirm that simulation is a necessary and important tool for constructing systems used in Satellite technology, particularly in the Attitude system of a CubeSat. The development of aerospace technology has a fundamental tool in artificial satellites and nanosatellites, including CubeSats. This same tool is available for countries with limited economies to carry out aerospace projects by building their CubeSats. This work is important because it provides the basis for future work and even the start of a project that can be sent into outer space.

Keywords

Magneto-Torquer, Control, PD, Attitude, CubeSat

1. Introduction

Nowadays, the use of CubeSat promotes the development of aerospace projects as their requirements facilitate their design and manufacturing for easy access to components and software, reducing the cost and the standard size that they must meet, which must be composed in 10 × 10 × 10 cms. Modules and a weight of no more than 1.33 kgs. The CubeSat cannot be compared with a conventional satellite in terms of performance and behavior as they came up with the necessity of designing, building, and operating smaller and more affordable satellites, letting research institutes in developing countries get involved in outer space.

A CubeSat has different controlling systems depending on the mission to be fulfilled, such as the Attitude, which positions it in outer space. There are three types of controlling systems of attitude that are used for its operation: chemical & electrical propulsion, reaction wheels, and magneto-torques. The Functional Principle of Magneto-Torquer consists of applying electric current to coils to carry a magnetic torque, which, along with the Ambiental magnetic field, produces a torque to rotate the CubeSat. This work focuses on the development of this type of Attitude system.

Considering the different types of control using inertial wheel systems, Magneto-torque, and other methods, there are studies on Attitude Control Systems. For example, as shown in [1], attitude control for CubeSat orientation is implemented through inertial wheels. The study describes the use of a microcontroller adapted to the characteristics of the system. To calculate the angular velocity and the direction of rotation of the wheels, four methods are used, applying mechanical simulations to model the movements of the CubeSat on a test platform. In [2], the design of a reaction or inertial wheel system is presented, where the angular momentum generated can be transferred in a single direction. This is possible thanks to the Möbius strip, which, after some modifications, becomes a surface capable of allowing the necessary reorientation, along with a secondary actuator that performs the displacement. With the simulation of the response of the motor elements and applying PID controllers, it is demonstrated that the design is effective, as it successfully achieves proper attitude regulation in finite time across all scenarios. In [3], two types of attitude systems were designed, fabricated, and characterized based on magnetic coils and inertial wheels for a 3U CubeSat. Together with the mathematical model of the satellite, the required torques for stabilization were determined, considering the limitations of mass, volume, power, and tolerance. In the attitude control simulation, the physical responses of the prototypes were verified to ensure they were correct, with small errors observed in both. In [4], the concept of exo-vela for a CubeSat is proposed using a parachute with an exobreak for maneuvers. PD controller feedback, with carefully adjusted gains, was proposed to prevent exceeding the control torques induced by the allowed exo-sail. The proposed concept was successfully demonstrated through simulations from the CubeSat control response, measuring the time it takes to correct the errors. The results of the attitude control proposed, and the strategy using the exo-sail are feasible.

To obtain information about the operation of the attitude control systems in the mentioned studies, simulations are conducted. However, none of these uses MATLAB-based software specifically designed for CubeSat simulations under outer space conditions, which in this study are considered constant and free from external or internal disturbances. For this reason, this study is considered both important and significant.

2. Literature Review

2.1. Description of Magneto-Torque and Operation, Its Theoretical Foundations

The Magneto-Torquer is a coil wound around with a conductor wire (usually copper) over some kind of support; this wire is coated with varnish to provide electrical insulation. This way of generating a magnetic field is known as an electromagnet, and the density generated due to the magnetic field may differ proportionally from the number of turns the coil might have and the quantity of current flowing through.

The coil is formed by fitting the core in the solenoid (Figure 1). Given that the iron has more permeability than the air, the lines of force can flow easily, and the magnetic field lines generated by the coil are concentrated in the core, producing a flow density increase inside the coil [5].

Figure 1. Representation of a current flowing through a solenoid and a coil and its lines of force.

2.2. Moment of Inertia a Solid

The moment of inertia of a solid is a scalar magnitude defined as [6]:

$I={\displaystyle \sum _{i=1}^N{m}_i{r}_i^2},$

$r_i$ represents the distance from the axis of rotation and $m_i$ as the masses.

Having said this, it can be deduced that the moment of inertia of a solid depends on the axis of rotation, and as mentioned before, when it is about a solid, the addition of the formula can be substituted by the following integral

$I=p{\displaystyle \int \text{d}m{r}^2},$

where, $\text{d}m$ represents the mass element of the solid directly related to the solid density $p$ and $r$ is the distance from the axis of rotation.

Given the mass element of the solid $\text{d}m$ and the solid density $p$ besides if this is homogeneous, it is possible to calculate the density of the integral when substituting $\text{d}m$ in the previous formula

$\text{d}m=p\text{d}V$

$I=p{\displaystyle \int r^2\text{d}V}$

where, $\text{d}V$ represents a solid volume and, to calculate a homogeneous solid’s inertia moment, it is necessary to solve the previous integral.

Taking the situation above into account, the variation per unit time is the linear velocity factor $\dot{r}=v$ and in case motion may take place along a radius circumference $r$ the magnitude of the linear velocity will be the following [7],

$v=r\frac{\text{d}\theta }{\text{d}t}=rw$ (1)

where, $v$ is the linear velocity, $w$ is the angular velocity.

Therefore, the angular moment $L$ applied to a rigid solid is obtained

$L={\displaystyle \int r\cdot mwr\text{d}m}={\displaystyle \int w{r}^2\text{d}m}=I\cdot w,$ (2)

it is known that a Magneto-Torquer operation is based on the ability to create a magnetic dipole; when a current flows through a coil, the equation of angular moment is obtained [8],

$L=I\cdot w.$ (3)

Remembering that the magnetic torque $\tau$ applied to a coil exposed to an external magnetic field may be expressed as the vector product of dipole moment and external magnetic field [5] [9],

$\tau =\mu \times B,$ (4)

where, $\mu$ is the dipole moment produced by the Magneto-torque, $B$ is the external magnetic field in such a way that the modulus of the magnetic torque can be expressed [9],

$\tau =\mu \cdot B\cdot sen\left(\theta \right),$ (5)

where $\theta$ is the angle that forms the axis of action of the magnetic force produced from the Magneto-Torquer and the desired angle. By saying this, it is possible to find the equation of motion of the system. To carry this on, it is necessary to apply the theorem of conservation of angular momentum. The law of conservation of angular momentum states if the net external torque is zero (it does not imply the external forces are zero, but it is an isolated system), the total angular momentum is conserved in other words, it is constant [6].

${\displaystyle \sum \text{Τ}}=I\cdot \dot{w},$ (6)

with $\text{Τ}$ , it is the torque applied to the system.

However, the magnetic torque produced by the magneto-torque, Equation (5), is not the only one that interacts with the system; there is also a resistive torque produced by different factors. From Equation (6), it can be deduced that when applying a torque to a solid, it can create an angular acceleration about the axis around which the object rotates, leading to the motion equation for the Magneto-Torquer

$\tau -{\text{Τ}}_r=I\cdot \omega ,$ (7)

where ${\text{Τ}}_r$ , represents the resistive torque which opposes the motion. Besides, the operating principle of a Magneto-Torquer is based on the torque produced by the magnetic dipole moment. Thus, it can possibly be related to the current intensity flowing through the coils along the satellite movement.

And so, relating the Equations (5), (7) it may bring

$\mu \cdot B\cdot sen\left(\theta \right)- {\text{Τ}}_r=I\cdot \omega ,$ (8)

the dipole moment produced by a coil is calculated as Equation (7) shows. The angular velocity function $\omega$ is obtained as:

$\omega =\frac{\mu B}{I}\cdot sen\left(\theta \right)-\frac{{\text{Τ}}_r}{I},$ (9)

where, $I$ represents the moment of inertia, $\mu$ is the dipole moment, $B$ is the external magnetic field, ${\text{Τ}}_r$ is a resistive torque, $\theta$ is the angle that forms the axis of action of the magnetic force produced by the Magneto-Torque, and the reference angle.

The inertia matrix is defined as a symmetric matrix whose elements are formed from the moment of inertia of three perpendicular axes and three inertia products. If $I\in {ℝ}^{3\times 3}$ respect to an origin $O$ , it may bring up [10] [11],

$I=\left[\begin{array}{ccc}{I}_x& -I_{xy}& -I_{xz}\\ -I_{xy}& I_y& -I_{yz}\\ -I_{xz}& I_{yz}& I_z\end{array}\right]$ (10)

where, $I=I^{\text{T}}$ .

Well, if the body’s coordinate system axes coincide with the main axes of inertia, that is, the principal frame of the satellite is chosen where the inertia matrix $I$ is diagonal, just like it happens with most of the satellites, the matrix is reduced as [10] [11],

$I=\left[\begin{array}{ccc}{I}_x& 0& 0\\ 0& I_y& 0\\ 0& 0& I_z\end{array}\right]$ . (11)

Therefore, the satellite’s Attitude dynamics are expressed by applying Euler’s equation:

$I_x{\omega }_x=-\left(I_z-I_y\right){\omega }_y{\omega }_z+{\tau }_x,$ (12)

$I_y{\omega }_y=-\left(I_x-I_z\right){\omega }_x{\omega }_z+{\tau }_y,$ (13)

$I_z{\omega }_z=-\left(I_y-I_x\right){\omega }_y{\omega }_x+{\tau }_z,$ (14)

where, $\omega =\left({\omega }_x,{\omega }_y,{\omega }_z\right)\in {ℝ}^3$ is the angular velocity vector of the satellite and $\tau =\left({\tau }_x,{\tau }_y,{\tau }_z\right)\in {ℝ}^3$ is the torque vector performing on the system.

3. Disruptions to be Considered

The CubeSats are normally designated for low Earth orbit missions, and as the orientation of the satellites of low orbit defines (Low Earth Orbiters, LEO), they are vehicles located in nearly circular orbits with altitudes between 300 km and 1500 km, which is an advantage due to the low altitude it is required a less power to reach their orbits compared to the conventional satellite. However, because of their orbit, these satellites may be affected by some perturbations due to the constant changes in atmospheric conditions. These perturbations usually affect the orientation and the orbit, which is why such disruptions must be considered in the Attitude control.

3.1. Geomagnetic Field

The geomagnetic field is the magnetic field extending from the Earth’s core. The Earth’s magnetic field may be approximated as a dipole inclined to 15 degrees to the axis of rotation of the Earth creates a magnetic field; its magnitude on the Earth’s surface has an average from 25 to 65 μT (microteslas) being stronger at the poles and weaker at the magnetic equator.

Table 1 shows the characteristics of the Earth’s geomagnetic field according to [12].

Table 1. Characteristics of the geomagnetic field.

Dipole field intensity	0.306 Gauss-Re3
Displacement of 3 dipoles	0.076 Re
Intensity of the Surface field (1 Re)	0.24 - 0.66 Gauss
Magnetic poles	Model WMM2020
Geocentric dipole	80.65 N, 72.68 O
Magnetic North pole	86.50 N, 164.04 E

where Re, defines the radius of the Earth’s model as 6.378 km. The dipole intensity is the dipolar portion intensity of the planetarium magnetic field outside the planet. In Gauss-Re³, Re is in units of the planet’s radius. (Dividing by the distance Re³ yields the field in Gauss).

3.2. Atmospheric Conditions

The Earth’s atmosphere is the gaseous layer of the Earth, being the most outermost and the least dense of the planet. This layer is composed of several gases that vary as the pressure at different altitudes. This mixture of gases that compose the atmosphere is genuinely called air. 75% of the atmospheric mass is at the first 11 km of altitude above the sea level. The main gases that compose it are Oxygen (21%) and nitrogen (78%), followed by argon, carbon dioxide, and water vapor.

Table 2 and Table 3 point out the characteristics and composition of the Earth’s atmosphere according to [12].

Table 2. Characteristics of the atmosphere.

Surface pressure	1014 mb
Surface density	1.217 kg/m2
Altitude scale	8.5 km
Atmosphere’s total mass	5.18 ×10+18
Average of temperature	288˚ k (15˚C)
Daily temperature rate	10˚C and 20˚C
Wind speed	0 and 100 m/seg
Average molecular weight	28.97

Table 3. Atmospheric composition (by fresh air volume).

Major gases	78.08% Nitrogen (N2), 20.9% Oxygen (O2).
Minor gases	Argon (Ar) 9340, Carbon dioxide (CO2) 410, Neon (Ne) 18.18, Helium (He) 5.24, Hydrogen chloride (HCV4) 1.7, Krypton (Kr) 1.14, Hydrogen (H2) 0.55.

4. Methodology

Below are some of the requirements and specifications a CubeSat model U1 must meet to be launched according to [13] the model mentioned above was selected for the simulations performed [14].

4.1. Mechanical Requirements

The CubeSats are cube-shaped satellites with the following specifications:

1) The coordinate system of a CubeSat 1U must coincide with the P-POD coordinate system (Poly Picosatellite Orbital Deployer), while integrated into the P-POD. The origin of the CubeSat system coordinates is located in the geometric center of the CubeSat.

2) The maximum mass of a CubeSat 1U is 1.33 kg.

3) The gravity center must be located no more than 2 cm. from the geometric origin in the X, Y & Z direction.

4.2. Electrical Requirements

1) The deployment switch must be operated all the time while it is integrated into the P-POD (Poly Picosatellite Orbital Deployer).

2) The CubeSat must have an RBF pin (Remove Before Flight).

3) The CubeSats will incorporate protection in the battery circuit to charge and discharge, thus avoiding unbalanced cell conditions.

5. PD Control

In this section, based on the position and angular velocity error a PD control system will be implemented.

Starting with the following equation

$u=K_p\theta a+K_d\left({\omega }_{orb}-{\omega }_{ob}\right),$ (15)

where $\theta$ and $a$ are Euler angles and axes, respectively, caused by a rotation given the matrix of attitude and the Euler axes. ${\omega }_{orb}$ is the orbital angular velocity with respect to the orbital system (angular velocity of the satellite). And ${\omega }_{ob}$ is the orbit angular velocity.

The PD control obtained is described as the following equation

$u=K_p\theta a+K_d{\omega }_{orb.}$ (16)

This control law is based on the Reijneveld and Choukroun algorithm, which is valid only if the satellite’s angular velocity is superior 4 times the revolutions per orbit. Otherwise, there will be absolutely no control [12].

Starting from

$M\left(q\right)\ddot{q}+C\left(q,\dot{q}\right)\dot{q}=\tau$ , (17)

where $q\in {ℝ}^n$ is the vector of generalized coordinates, $M\left(q\right)\in {ℝ}^{n\times n}$ is the inertia matrix, $C\left(q,\dot{q}\right)\dot{q}\in {ℝ}^n$ is the angular velocity vector, and $\tau \in {ℝ}^n$ is the torque vector acting on the system.

It will use the structure of a PD defined as follows

$\tau =K_p\tilde{q}-K_v\dot{q},$ (18)

$\tilde{q}=q-q_d\in {ℝ}^n$ is the position error, $q_d\in {ℝ}^n$ is the desired position, $K_p,K_v\in {ℝ}^{n\times n}$ are the proportional and derivative gains.

The state variables (canonical form of the controller) obtained when substituting the PD control into the dynamic model are

Position error:

$x_1=\tilde{q}=q-q_d.$ (19)

Note: When performing position control (with regulation), the desired position $q_d\in {ℝ}^n$ is assumed to be constant.

Velocity error:

${\dot{x}}_1=\dot{q}=x_2.$ (20)

Acceleration error

${\dot{x}}_2=\ddot{q}=M{\left(q\right)}^{-1}\left[K_p{x}_1-K_v{x}_2-C\left(q,\dot{q}\right)x_2\right]$ . (21)

Afterwards, the equation is derived as follows

$\left[\begin{array}{c}{\dot{x}}_1\\ {\dot{x}}_2\end{array}\right]=\left[\begin{array}{c}{x}_2\\ M{\left(q\right)}^{-1}\left[K_p{x}_1-K_v{x}_2-C\left(q,\dot{q}\right)x_2\right]\end{array}\right]$ (22)

$\left[\begin{array}{c}{\dot{x}}_1\\ {\dot{x}}_2\end{array}\right]=\left[\begin{array}{cc}0& I\\ M{\left(q\right)}^{-1}{K}_p& -M{\left(q\right)}^{-1}\left[K_v+C\left(q,\dot{q}\right)\right]\end{array}\right]\left[\begin{array}{c}{x}_1\\ x_2\end{array}\right].$

Once the equation homogenous ${\dot{x}}_1=x_2=0$ , the existence, uniqueness, and location of the equilibrium point are found out such that, $x_1=x_2=0$ , implying the equilibrium point is $\left(x_1,x_2\right)=\left(0,0\right)$ ,

$\left[\begin{array}{c}0\\ 0\end{array}\right]=\left[\begin{array}{c}{x}_2\\ M{\left(q\right)}^{-1}{K}_p{x}_1-M{\left(q\right)}^{-1}\left[K_v+C\left(q,\dot{q}\right)\right]x_2\end{array}\right]$ (23)

Stability

Stability can be tested with the Lyapunov function criterion [15]. Considering the following positive define function (Lyapunov candidate function),

$V\left(x_1,x_2\right)=\frac{1}{2}{x}_1^{\text{T}}{K}_p{x}_1+\frac{1}{2}{x}_2^{\text{T}}M\left(q\right)x_2,$ (24)

the representation in terms of position and velocity errors is obtained by substituting the state variables in Equation (24), and the result is a positively defined function

$V\left(\tilde{q},\dot{q}\right)=\frac{1}{2}{\tilde{q}}^{\text{T}}{K}_p\tilde{q}+\frac{1}{2}{\dot{q}}^{\text{T}}M\left(q\right)\dot{q},$ (25)

when deriving Equation (25):

$\dot{V}\left(\tilde{q},\dot{q}\right)={\tilde{q}}^{\text{T}}{K}_q\dot{\tilde{q}}+{\tilde{q}}^{\text{T}}M\left(q\right)\dot{\tilde{q}}+\frac{1}{2}{\dot{q}}^{\text{T}}\dot{M}\left(q\right)\dot{q},$ (26)

The state variable $x_2$ in Equation (21) will give:

$M\left(q\right)\ddot{q}=K_p\tilde{q}-K_v\dot{q}-C\left(q,\dot{q}\right)\ddot{q},$ (27)

when substituting Equations (26), (27):

$\dot{V}\left(\tilde{q},\dot{q}\right)={\tilde{q}}^{\text{T}}{K}_p\dot{\tilde{q}}+{\dot{q}}^{\text{T}}{K}_p\tilde{q}-{\dot{q}}^{\text{T}}{K}_v\dot{q}+\frac{1}{2}{\dot{q}}^{\text{T}}\left[\dot{M}\left(q\right)-2C\left(q,\dot{q}\right)\right]\dot{q}$. (28)

By applying the antisymmetric matrix properties in Equation (28), we get:

$\dot{V}\left(\tilde{q},\dot{q}\right)={\tilde{q}}^{\text{T}}{K}_p\dot{\tilde{q}}+{\dot{q}}^{\text{T}}{K}_p\tilde{q}-{\dot{q}}^{\text{T}}{K}_v\dot{q}.$ (29)

When using the state variable ${\dot{x}}_1$ , Equation (19) will give:

${\dot{x}}_1=\dot{\tilde{q}}=\{\begin{array}{l}\dot{q}-{\dot{q}}_d=\dot{q}\\ {\dot{q}}_d-\dot{q}=-\dot{q}\end{array}$ (30)

And applying the properties of the transpose

${\left(xyz\right)}^{\text{T}}=x^{\text{T}}{y}^{\text{T}}{z}^{\text{T}},$ (31)

${\left(x^{\text{T}}\right)}^{\text{T}}=x,$ (32)

considering that $K_p$ is symmetric, as, $K_p=K_p^{\text{T}}$ and with the described properties in Equations (31), (32), (30) can be modified as follows

$\dot{V}\left(\tilde{q},\dot{q}\right)=-{\dot{q}}^{\text{T}}{K}_v\dot{q}.$ (33)

Finally, the derivative of the function of Lyapunov is less or equal to zero; this is because when the error of the position is zero, this term $-{\dot{q}}^{\text{T}}{K}_v\dot{q}$ becomes zero. In other words, it is indeed less or equal to zero. Concluding that the equilibrium point $x_1=x_2=x^{*}\left(0,0\right)$ located at the origin of the state space, is stable.

6. Simulation and Results

The Propat functions can be grouped into transformations of coordinates, time and ephemerides, transformations and propagation of Attitude, and environmental property conversions and orbit propagation. The sensors used for the Propat are magnetometers of 3 axes and coarse solar sensors along 3 orthogonal magnetic torsion coils set for Attitude control [16]. The main coordinated systems used by the Propat are the geocentric inertial system, the geocentric terrestrial system, and the orbital system. The satellites use sensors to define their location in space and their orientation. These sensors work based on the properties of the outer space environment, such as the Earth, Sun, Moon, stars, and the geomagnetic field, among others directions.

This is why it is important to always simulate these sensors, the Propat library currently includes 3 models: the Sun’s position in the geocentric inertial system, the Earth’s orientation at a given time, and the Earth’s magnetic field. Additionally, there is a function to demonstrate if the satellite is in the Earth’s shadow.

6.1. Variables for Simulation

For this case, the face X of the satellite is oriented in order to calculate the Attitude matrix between the satellite’s frame and the orbital system based on the Attitude vectors function of the Sun and the geomagnetic field measured by the sensors. In addition, a specific orbit over the State of Puebla was chosen.

Three sets of orthogonal magnetic torsion coils are used for the Attitude control, and for this case, it is used a NCTR-M002 Magneto-Torquer (see Table 4) whose characteristics make it particularly useful to be implemented for the CubeSat due to its low use of energy and its lightweight. Table 5 shows the simulation parameters used.

It was selected as a 1U CubeSat model as shown in Figure 2.

Table 4. NCTR-M002 magneto-torquer rod.

Magnetic moment	>0.24 m2
Functionality Range	−35˚C and 75˚C
Power supply voltage	5 V
Mass	<30 grs

Table 5. Simulation parameters.

Initial Euler angles	[30 50 10]
Initial angular velocity	[0.1 0 0.5] rad/seg
Orbit	Polar (90˚)
Eccentricity	0.0002
Altitude	300 km
Argument of perigee	111.196˚
Average anomaly	248.196˚
Ascending node longitude	6.524 hrs

Figure 2. Model of the 1U CubeSat selected.

With a mass equal to M = 1.2 kg and with an inertia matrix

$I=\left[\begin{array}{ccc}2.166\text{e}-3& 0& 0\\ 0& 2.166\text{e}-3& 0\\ 0& 0& 2.166\text{e}.3\end{array}\right]$ (34)

6.2. Elements for the Simulation

The next Keplerian elements were taken into consideration for the simulation, such as the Earth’s radius, orbit, eccentricity, and the orbit’s inclination. For this case, a polar orbit was selected, and the rest of the elements, such as the perigee argument, average anomaly, and ascending node longitude, were selected.

The variation of the Keplerian elements is calculated according to the Earth’s shape. An initial Attitude (initial position) for the CubeSat in radians is proposed, and an angular velocity is given in radians per second.

$eulxyz={\left[30,50,10\right]}^{\prime }\frac{\pi }{180};$

$w_{ang}={\left[0.1,0,0.5\right]}^{\prime };$

The inertia matrix is symmetric as it is with most of the satellites, the coordinate system axes coincide with the principal inertia axes. The matrix values must be given in Kg.

Three magnetic torsion coils are used for the Attitude control, one for each axis, with a magnetic moment of 0.2 Am².

The satellites usually use sensors to find out their location and orientation in orbit (Attitude). These sensors work based on the direction of the Earth, Sun, Moon, stars, the geomagnetic field, and among others. This is why it is important to keep the sensors mentioned during the simulation in mind. Propat offers distinct models that allow simulating the Sun’s position in the geocentric inertial system, the Earth’s orientation in a specific moment (Greenwich Sidereal Time), and the Earth’s magnetic field. The simulation code shows exactly the model used to test the sensors at a specific time.

The time propagation of the simulation is given in seconds. The time, angles, angular velocity, and the orbit vector start initializing. The control torque and the current of the magneto torques start initializing.

As a numeric integration is estimated, the precision is defined. With this, it is now possible to calculate the Orbit and Attitude propagation. Also, the propagation of the orbit and conversion of the Keplerian elements to a vector of state.

In this case, perturbations will not be considered. Therefore, the Attitude vector will be defined and initialized.

6.3. X Face of the Satellite Oriented towards the Earth

It is necessary to have knowledge of the direction and the intensity of the geomagnetic field to orient a Satellite’s face towards the Earth.

It is also recommended that the ODE45 numerical integrator be used for the propagation of the Attitude.

Additionally, it is necessary to know not only the position vector but also the Sun’s position vector in the Satellite’s frame.

The Gaussian was added to the Sun’s position vector and to the geomagnetic field since the data of the sensors is corrupted by the line noise before being processed, thus limiting the transmission distance.

The axes are adjusted so that the X axis aims at the Earth, the Y axis is orthogonal to the X axis in the orbital plane, and the Z axis is normal to the orbital plane.

Considering two reference systems, one belonging to the Earth and the other belonging to the satellite, the Attitude matrix is calculated for two pairs of vectors known in two different reference systems. If the two coordinate systems are the same, the matrix gets its identity back, as shown in Figure 3. At last, the orbital angular velocity is calculated. As was mentioned previously, the control is valid only if the satellite’s angular velocity is a minimum of four times the revolutions per orbit.

Figure 3. Coordinate system.

6.4. The Control Code in Use during the Simulation

The current supplied to the Magneto-Torquer is calculated considering the number of spins and the maximum current they can hold.

Current in the coils

EB = teta’.*earth_field;

K = 1e5;

MI = −k*(cross(EB, (w_ang-wo_b)));

N = I max;

xyzC = MI/(710*0.04);

xcurrent = xyzC(1);

ycurrent = xyzC(2);

zcurrent = xyzC(3);

Finally, the data is saved for plotting.

time = cat(2, time, t);

omeg = cat(2, omeg, w_ang);

euler = cat(2, euler, teta*180/pi);

corx = cat(2, corx, xcurrent);

cory = cat(2, cory, ycurrent);

corz = cat(2, corz, zcurrent).

6.4.1. Simulation Results of the PD Control and Discussions

According to the next expression of the control PD, we can say:

$u=K_p\theta a+K_d\left({\omega }_{orb}-{\omega }_{ob}\right),$ (35)

$\theta$ and $a$ are Euler’s angle and axis, while ${\omega }_{orb}$ is the orbital angular velocity with respect to the orbital system and ${\omega }_{ob}$ is the orbital angular velocity. Meanwhile, the selection of the gain $K_p=-0.00025$ and $K_d=-0.002$ was done based on the criteria of $K_p<K_d$ since this is the way to prioritize the equilibrium point, however, if the value of the gains is reversed, the system reduces its angular velocity, making the system reach its equilibrium point in a longer period [16].

% Vangular-Vorbital

cont_torq = −0.00025*teta − 0.002*(w_ang-wo_b)

6.4.2. Simulation Graphics

Considering the gains of $K_p$ & $K_d$ and the characteristics of the satellite the following results are obtained, which represent the orientation measured in Euler angles, the angular velocity measured in rad/s, and finally, the variation of the current supplied to the Magneto-Torquer.

As Figure 4 shows, it can be observed that the satellite reached the desired position even after the position error was zero due to the line noise.

Figure 4. Position error.

In Figure 5, it can be clearly observed that when the satellite reaches the desired position, the angular velocity is equal to zero.

Figure 5. Angular velocity.

As it is shown in Figure 6, the current supplied does not exceed 0.03 A because the maximum current the Magneto-Torquer can handle is 0.04 A.

Figure 6. Current supplied to the Magneto-Torque.

7. Conclusions

The simulation carried out for the Attitude control of a CubeSat, using the control PD, provides significant information for the already mentioned mini-satellites development and the implementation of the control system. As we know, the development of this technology is based on high budgets, which, unfortunately, we may not be able to afford due to different circumstances. The simulation is a tool that does not require a high budget to obtain the study and analysis results.

The results that were obtained show us precisely the data gathered considering the conditions of the outer space the CubeSat was exposed to, the variables that must be kept in mind, and its consideration for the development and implementation, which could be sent to the outer space.

For the development of this knowledge applied to Satellite Technology, it is important to consider the national context, whether we are slightly or not supported, and notice the difference that exists in many countries for these applications and their importance for their development in every country. Another reason to consider the importance of the work is to recognize how many future projects will continue to keep this line of investigation in demand.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Lechuga, V.M. (2017) Implementación del Control de Actitud de un CubeSat mediante Ruedas de Reacción. https://www.researchgate.net/publication/359045661_Implementacion_del_Control_de_Actitud_de_un_CubeSat_mediante_Ruedas_de_R[eaccion
[2]	Lechuga, V.M. (2023) Diseño de un Arreglo Reconfigurable de Ruedas de Reacción basado en la Banda de Möbius para el Control de Actitud de Satélites CubeSat. https://www.researchgate.net/publication/370952991_Diseno_de_un_Arreglo_Reconfigurable_de_Ruedas_de_Reaccion_basado_en_la_Banda_de_Mobius_para_el_Control_de_Actitud_de_Satelites_CubeSat
[3]	Prado-Morales, J., Navarrete-Montesinos, M., Molina, J.P. and Hernández-Arias, H. (2019) Diseño e implementación de actuadores del subsistema de determinación y control de orientación de un satélite CubeSat 3U. https://www.researchgate.net/publication/372855243_Diseno_e_implementacion_de_actuadores_del_subsistema_de_determinacion_y_control_de_orientacion_de_un_satelite_CubeSat_3U
[4]	Shan, S. and Westfall, A. (2019) Attitude Control System Design for CubeSats Configured with Exo-Brake Para Chute. https://ntrs.nasa.gov/api/citations/20190001361/downloads/20190001361.pdf
[5]	Hermosa, D.A. (1999) Principios de Electricidad y Electrónica II. Marcombo, 50-62. https://www.academia.edu/41549740/Principios_de_Electricidad_y_Electronica_I_Antonio_Hermosa_Donate
[6]	Landau, L.D. and Lifshitz, E.M. (1970) Física teórica. Mecánica. Reverte, 15-28. https://www.reverte.com/libro/fisica-teorica-mecanica-vol-i_91625/
[7]	Marion, J.B. (2014) Dinámica clásica de las partículas y sistemas. Reverte, 409-455. https://books.google.co.ve/books?id=MpblAXz_StkC&printsec=frontcover#v=onepage&q&f=false
[8]	Tipler, P.A. and Mosca, G. (2012) Física para la ciencia y la tecnología. Volumen 1A Mecánica. Reverte, 331-354. https://www.academia.edu/45147942/F%C3%ADsica_Tipler_6a_Ed_Vol_1
[9]	Tipler, P.A. and Mosca, G. (2004) Física para la ciencia y la tecnología. Volumen 2 de Electricidad y magnetismo/Luz. Reverte, 887-903. https://archive.org/details/fisica-para-la-ciencia-y-la-tecnologia-by-paul-a.-tipler-gene-mosca-vol-2-electr
[10]	Sidi, M.J. (1997) Spacecraft Dynamics and Control: A Practical Engineering Approach. Cambridge University Press, 22-92. https://www.google.com/search?q=Sidi%2C+M.+J.+(1997).+Spacecraft+Dynamics+and+Control%3A+A+Practical+Engineering+Approach.&oq=Sidi%2C+M.+J.+(1997).+Spacecraft+Dynamics+and+Control%3A+A+Practical+Engineering+Approach.+&gs_lcrp=EgZjaHJvbWUyBggAEEUYOdIBCTI2ODZqMGoxNagCCLACAQ&sourceid=chrome&ie=UTF-8
[11]	Desouky, M.A.A., Prabhu, K. and Abdelkhalik, O. (2018) On Spacecraft Magnetic Attitude Control. 2018 Space Flight Mechanics Meeting, Kissimmee, 8-12 January 2018, 4-6. https://arc.aiaa.org/doi/10.2514/6.2018-0205 https://doi.org/10.2514/6.2018-0205
[12]	Williams, D.R. (2020) Earth Fact Sheet. 15/03/2021, de NASA Goddard Space Flight Center. https://nssdc.gsfc.nasa.gov/planetary/factsheet/earthfact.html
[13]	NewSpace (2016) Varilla Magnetorquer NCTR-M002. 12/01/2021, de CubeSatShop, 7-44. https://www.cubesatshop.com/product/nctr-m002-magnetorquer-rod/
[14]	Munakata, R., Carnahan, J. and Pignatelli, D. (2014) CubeSat Design Specification. California Polytechnic State University, 9-15. https://static1.squarespace.com/static/5418c831e4b0fa4ecac1bacd/t/62193b7fc9e72e0053f00910/1645820809779/CDS+REV14_1+2022-02-09.pdf
[15]	Ogata, K. (2010) Ingeniería de control moderna. Pearson Education. https://www.academia.edu/9814191/Ingenieria_de_Control_Moderna_Ogata_5ed
[16]	Carrara, V. (2015) An Open Source Satellite Attitude and Orbit Simulator Toolbox for Matlab. Proceedings of the XVII International Symposium on Dynamic Problems of Mechanics, Natal, 22-27 February 2015, 2-8. https://www.researchgate.net/publication/274065310_An_Open_Source_Satellite_Attitude_and_Orbit_Simulator_Toolbox_for_Matlab