Vision-Based CubeSat Closed-Loop Formation Control in Close Proximities

4.1 Pose estimation methodology

Figure 2 shows the leader and follower CubeSats. Two poses are shown for the follower. The one on the left is the desired (ideal) position of the follower CubeSat with respect to the leader CubeSat. However, the follower has drifted away from its desired pose with respect to the leader. The actual pose of the follower after the drift (shown on the right) is denoted by a spatial frame {f}.

Fig. 2

The leader and follower CubeSats with body frames

Note that the transformation matrix Trfd is known from Eq. (1). If one can determine the transformation between the frames {f} and {f^d} (i.e., the matrix Tfdf), then, the actual relative pose of the follower {f} with respect to the leader {r} can be determined by the following simple transformation.

The Mirage pose estimation method directly calculates the elements of matrix Tfdf via an analytical linear non-iterative solution, when two (or more) cameras are used.

The calculation of Tfdf involves the comparison of two images: the desired image and the actual image, which are defined in the following. Four non-planar markers (with known geometry) are installed on the leader (points 1 through 4 in Fig. 2). Two cameras are installed on the follower, but only one of them, represented by frame {c_i}, is shown in Fig. 2. The camera {c_i} (i = 1, 2) generates images of the markers. The image seen by the camera when the follower is at its desired pose is called the desired image. On the other hand, the image that the camera sees when the follower is in its actual pose is called the actual image.

The difference in the pixel coordinates of marker j (j = 1, …, 4) in the actual and desired images of camera c_i (i = 1, 2) is a function of Tfdf. This difference in pixel coordinates is denoted by the following equation.

where (x_pj, y_pj)_ci and (xpjd,ypjd)ci are the pixel coordinates of marker j in camera i in the actual and desired image, respectively. Note that both the actual and the desired pixel coordinates are known. The actual pixel coordinates are found by extracting the pixel coordinates of the markers in the actual image via image processing. The desired pixel coordinates are found as follows.

After markers are physically installed on the leader, the Euclidean position vector of marker j with respect to the origin of {r} is denoted by (r_r)_j = [rx,ry,rz,1]jT. The desired 3D projection of marker j in the camera space (pcid)j=pcixdpciydpcizdjT is found by a series of transformation applied to the vector r_r.

Here, the 4 × 4 matrix Tfdci transforms any vector from frame {f^d} to frame {c_i}, and the 3 × 4 matrix K_ci contains the intrinsic camera parameters (focal length, principle points and distortion parameters), which are available a priori via camera calibration. Then, the desired pixel coordinates are obtained.

Therefore, the pixel error ejci is known. If a relation between the pixel error ejci and the transformation matrix Tfdf is established, that relation can be used to solve for the elements of Tfdf. The derivation of such a relation is explained in the following.

First, the pixel error ejci can be calculated based on the coordinates in the 3D camera space.

The actual 3D projection of marker j in the camera space {c_i} is denoted.

Equation (7) is rewritten by incorporating the fact that Trf=TfdfTrfd.

Note that there are no unknowns in Eq. (4) and the only unknowns in Eq. (8) are the twelve elements of the matrix Tfdf.

Equation (9) is substituted in Eq. (8). The result of the substitution and Eq. (4) are expanded to obtain the components of (p_j)_ci and (pjd)ci. The resulting components are substituted in Eq. (6). After very involved manipulations, which can be found in [15], the following relation is found between the pixel error ejci and the elements of Tfdf.

where V^c_i is a 2 × 12 matrix and Y^c_i is a 2 × 1 vector with known values. For the detailed formulation of these parameters, see [15]. The 12 × 1 vector tfdf contains all 12 unknown elements of the transformation matrix Tfdf.

Note that Eq. (10) includes 2 equations and 12 unknowns and cannot be solved in its current form. However, this equation is valid for camera c_i (i = 1, 2) and marker j (j = 1, …, 4). For this case, Eq. (10) is written 8 times, covering the combinations of 2 cameras c₁ and c₂ and 4 markers.

Equation (12) forms a linear system of 16 equations in terms of 12 unknown elements of tfdf, which can be easily solved in the Least-Square sense.

Now that the vector tfdf is calculated, the transformation matrix Tfdf is found. Then, Tfr is calculated using Eq. (2). The elements of Tfr are related to the pose parameters of the follower (q = [x, y, z]^T and Φ = [ϕ, θ, ψ]^T).

Therefore, knowing Tfr one can solve for q = [x, y, z]^T and Φ = [ϕ, θ, ψ]^T.

4.2 The simulated cameras

The block “simulated cameras” calculates the “actual pixel coordinates” of the marker j(j = 1, …, 4) for simulation only using the simulated pose of the follower CubeSat mathematical model. The actual 3D projection of marker j in the camera space is found using Eq. (7) and Trf=(Tfr)−1. Finally, the simulated pixel coordinates are obtained as follows.

6.1 Governing equations for simulations

The coordinate axes shown in Fig. 2 are considered. Let xrr and xfr denote the position of the leader and the follower satellites with respect to the geocentric system (X, Y, Z), but expressed along the moving axes (x, y, z) of frame {r} at the location of the leader satellite. Note that axes x, y, and z are radial, along track, and cross-track directions, respectively, defined by Hill and used in [2].

It is known that the effect of the follower’s thrusters force u_f on the acceleration of the follower x¨fr is described by the follower’s equation of the motion.

In this equation, m_f is the mass of the follower satellite; FEfr is the pulling force from the Earth on the follower satellite; Frfr is the pulling force from the leader satellite on the follower satellite; and Fdfr is the resultant perturbing forces on the follower satellite. All vectors, except u_f, have a superscript r, which indicates that they are written along the axes of the moving coordinate (x, y, z) of the frame {r}. The thrust vector u_f is expressed in frame {f}, and Rfr transforms a vector from frame {f} to frame {r}.

This equation is integrated for simulations in the “formation dynamic model” block to provide xfr. The trajectory of the leader xrr is directly defined by the user.

6.2 Governing equations for control law derivation

Here, to achieve a precise formation, the relative position of the follower with respect to the leader will be controlled. Therefore, the dynamic equations of motion must be written as a relation between q̈ and the control forces from the follower’s thrusters u_f.

The acceleration of the follower satellite can be written in terms of the acceleration of the leader satellite.

where q is the relative distance of the follower with respect to the leader along the moving axes (x, y, z) of frame {r}, and ωrr is the instantaneous rate of rotation of the leader around the Earth expressed along the axes of {r}. For the derivation of the controller, the “nominal” situation is used, in which the resultant perturbing forces on the follower satellite are neglected (Fdfr=0) and the nominal mass of the satellite is used (m̂_f).

Equation (22) is substituted in Eq. (23), and the result is solved for q̈. The following is obtained

where

in which F^Efr and F^rfr are calculated using m̂_f.

Equation (24) describes the direct effect of the follower’s thrust u_f on the relative motion of the follower with respect to the leader (q). This equation will be used for derivation of the formation controller.

6.3 Implementation notes

A controller that is derived based on Eq. (24) needs the measurement of the realtime acceleration of the leader (x¨rr), the relative position of the follower with respect to the leader (q), and its rate (q̇). The acceleration x¨rr can be measured by an accelerometer on board of the leader satellite, and can be broadcast to the follower satellite. The relative position q is measured by the vision system, and its rate q̇ is observed by the HOSM differentiator.

7.1 Governing equations of motion

The ability to correct the satellite’s attitude is important, because after long periods of time, the external perturbing forces and moments will eventually rotate the satellite bodies. This may cause the loss of sight for vision-based formation correction maneuvers. Therefore, an attitude controller must be used to realign the satellites precisely when necessary.

A CubeSat body axes (x_B, y_B, z_B) and its three orthogonal reaction wheels are shown in Fig. 3. The moment components τ = [τ_1x, τ_2y, τ_3z]^T are exerted by the motors of the reaction wheels on the CubeSat’s body. The following equations describes the rotational dynamics of the CubeSat body.

Fig. 3

A CubeSat with three orthogonal reaction wheels

The rotational dynamics of the reaction wheel are described by the following equation, in which θ_i(i = 1, 2, 3) is the angular position of reaction wheel i.

The derivation and the definition of the notation used in the above equations are found in the Appendix.

In summary, Eqs. (27), (28), and (29) are the complete dynamic model of the CubeSat and reaction wheels as a system, which are inserted in the “attitude dynamic model” block in Fig. 1. At any simulation step, the reaction wheel driving moments are known. One can calculate the CubeSat’s angular acceleration vector ω̇_B using Eq. (27). Also, the reaction wheel driving moments and the CubeSat’s angular acceleration vector ω̇_B just obtained from Eq. (27) can be used with Eq. (29) to calculate the reaction wheel’s angular acceleration (θ̈₁, θ̈₂, θ̈₃). Integrating these angular accelerations will provide ω_B and (θ̇₁, θ̇₂, θ̇₃) for the next time step. The integration can go on by repeating the above steps.

7.2 Governing equations for control law derivations

Equation (27) is written in a more concise form.

where

are calculated using the “nominal” values of the system parameters. The following outputs are defined for the control system in Eq. (30).

where Φ^d is the desired orientation of the satellite, which in general can be a function of time. A controller must be derived that can stabilize the output y of the system in Eq. (30) and its time derivative at zero vector in finite time. The dynamics of the output y is derived for control design in the following.

The second time derivative of the surface variable is calculated.

If the notations f̂_ϕ = Ef̂_ω + Ėω_B – Φ̈^d and b̂_ϕ = Eb̂_ω are adopted, the model for attitude control derivation becomes.

The above equation is used for the derivation of a HOSM attitude controller.

8.1 Nominal control

A change of state variables such as σ = z₁ and σ̇ = z₂ and change of the control signal such as ŵ = f̂ + b̂û are applied to the output system in Eq. (35). The following chain form is obtained.

It can be shown that the following feedback control law can stabilize the system in Eq. (36) at zero vector in finite time [16],

where v_1j = 3/5 and v_2j = 3/4, and the control gains k_ij are found such that the following characteristic polynomial (p² + k_2jp + k_1j = 0, j = 1, 2, 3) is Hurwitz. The above formulation is used as explained in the following. At each time step, the measurement of the real-time errors and their rates provide the feedback signals z₁ and z₂.

Then, the nominal control is calculated using Eq. (37). Finally the equivalent control is calculated by inverting the relation ŵ = f̂ + b̂û.

The above “nominal” control works accurately if the system model has no uncertainty, that is if f̂ and b̂ are exact. However, normally, the real system is represented by {f} and {b}, which are different due to modeling and parameters uncertainties. In the following, the nominal control is amended to account for such uncertainties.

8.2 Extension to uncertain systems

The relation between the terms of the nominal system model and those of the actual system are defined as f = f̂ + f̃ and b = b̂ + b̃, in which f̃ and b̃ are bounded. The following equation represents the real system.

Noting that z₂ = σ̇ and substituting Eq. (40) in the above relation, one obtains.

To make sure that ż₂ → 0 even in the presence of uncertainty, the nominal control ŵ in Eq. (42) is replaced by

where s = z₂ + z_aux, and G is a positive constant. The gain G must be determined such that it is large enough to compensate for uncertain terms. This gain is determined using a Lyponuv design method.

8.3 Lyaponuv gain design

A Lyaponuv function is defined as V = 12s^Ts. The gain G must be determined such that V̇ = s^Tṡ ≤ 0 for all s. Because ṡ = ż₂ + ż_aux = ż₂ – ŵ, V̇ is derived as follows.

To design the controller gain G, the following bounds are defined for the uncertain terms.

Note that α ≈ 1. Knowing that c^Td ≤ ∥c∥ ∥d∥ and ∥s∥ ≤ s^Tsign(s), and using Eq. (45), one can obtain the following.

Equation (47) implies that, if the control gain G is selected to meet the following requirement,

where η > 0, then, V̇ ≤ –η ∥s∥, or V̇ is negative for all s and V = 12s^Ts vanishes to zero despite bounded uncertainties in the system defined by Eq. (46), provided that the following control law is used.

Based on the formation dynamic model in Eq. (24), the above control law is used as the “HOSM formation controller” by setting f̂ = f̂_x, b̂ = b̂_x, z₁ = q, z₂ = q̇, and u = u_f. Also, based on the attitude dynamic model in Eq. (34), the above control law is used as the “HOSM attitude controller” by setting f̂ = f̂_ϕ, b̂ = b̂_ϕ, z₁ = y, z₂ = ẏ, and u = τ. The term ŵ is always calculated using the generic Eq. (37).