Printer Friendly

Consideration of Three-Dimensional Attitude and Position Control for a Free-Floating Rigid Body Using Three Thrusters.

1. Introduction

What is the minimum number of thrusters required to control both the attitude and position of a rigid satellite? What strategies can be applied and what kind of motion is controlled for a satellite in which there are few remaining thrusters? These questions motivated us to pursue the current investigation because nonlinear underactuated control makes these satellites fault-tolerant with respect to thruster malfunctions and enables continuity of the respective missions.

Most satellites are equipped with a sufficient number of thrusters to control the attitude and position in three dimensions because of the following reasons. In typical thrusters, the force directions are fixed to a satellite body and the magnitudes are restricted as positive. This implies that another thruster is required to generate a force in the opposite direction. Consequently, more than twelve thrusters are required for a satellite to guarantee the six degrees of freedom in motion and to ensure the robustness of the system in simple control logic. However, the number of thrusters can be reduced by applying nonlinear control methods. Hence, combined attitude and position control based on the nonlinear motion of a satellite can achieve a significant reduction in the number of actuators.

A few underactuated systems in which the degrees of freedom of motion exceed the number of actuators can be potentially completely controlled by utilizing their nonholonomic constraints [1, 2]. A free-floating satellite is a type of an underactuated system, and its attitude and position dynamics are represented as nonintegrable second-order differential equations that act as nonholonomic constraints [3]. Several intensive studies focused on the attitude motion control of a satellite by utilizing nonholonomic constraints. Several extant studies indicate that two independent torques successfully control a three-dimensional satellite attitude. Tsiotras and Longuski [4] proposed new attitude parameters termed as "wz-parameters," and then Tsiotras and Doumtchenko [5] developed a discontinuous control law that utilized the parameters. In another study, Morin and Samson [6] designed a smooth time-varying feedback control. Recent studies shifted to more practical ones including a controller for the underactuated CMG satellite system [7], implementation of the kinematic planning scheme [8], and an optimal process for state stabilization [9]. The strategies indicate that two bidirectional control torques (i.e., both positive and negative directional torques are generated) with independent magnitude can control three attitude parameters of a satellite. Therefore, the problem of satellite attitude control is transformed into the number of thrusters that are required to generate independent torques with respect to the two principal axes. Sidi [10] presented a pioneering study related to this problem that numerically demonstrated that four thrusters produce independent torques with respect to three principal axes. In [11], we analytically confirmed the forementioned results and indicated that three thrusters that satisfy a configuration produce bidirectional torques with respect to two independent axes. This implies that when the nonholonomic dynamics of a satellite are considered, three thrusters control the three-dimensional attitude of a satellite. Hence, Matsuno et al. [11] and Yoshimura et al. [12] demonstrated an attitude control strategy and developed feedback controllers although a zero or unidirectional (only positive or negative directional) torque was generated under the thruster configuration around the third axis owing to the coupling effect with the other two axes.

The present study focuses on the simultaneous attitude and position control of an underactuated satellite. In two-dimensional motion, at least two thrusters are required to accelerate and decelerate the rotation of a satellite. Essentially, the same type of controllers in [13-15] is applied to obtain a controller with two thrusters for the in-plane attitude and position motion of a free-floating satellite. The study numerically verified that a switching controller that utilizes an invariant manifold achieves a two-dimensional attitude and position control by using positive albeit variable magnitudes of thruster forces. The study was extended in [16] for two on/off types of thrusters in which the input magnitudes are fixed.

In contrast to the two-dimensional cases, only a few studies discuss three-dimensional attitude and position control for underactuated satellites. With respect to a satellite with four thrusters, the previous study [17] described a design procedure of motion planning for the simultaneous control of the attitude and position. However, in the case with three thrusters, the three-dimensional attitude and position control is quite challenging because of the following reason. Three thruster configurations that guarantee bidirectional control torques for two principal axes generate dependent torque with respect to the third principal axis of a satellite. Specifically, the attitude control proposed in the present requires a series of complex nonlinear processes. Furthermore, thruster forces are unidirectional, and thus, the controllability of the system is uncertain. It should be noted that although the theorem of Sussman [18] or Goodwine and Burdick [19] provides sufficient conditions for the controllability of nonholonomic systems, they cannot be applied in this study owing to the nonnegative restriction of the thruster forces.

Given the above background, the present study develops a switching controller by using three invariant manifolds and provides a motion planning scheme in the final phase. This paper deals with the problem to control a free-floating rigid body in three-dimensional space using three thruster forces (hereafter, we refer to this type of rigid body model as "FFRB"). In this simplified problem, the proposed method achieves both attitude and position control by using only three thrusters in which the directions are fixed with respect to the body. The effectiveness of the method is numerically certified in ideal conditions (constant mass and nonlimited positive thruster force). This study is the first to demonstrate that the three-dimensional attitude and position of a rigid body can be controlled with three fixed thrusters, although the following problems remain in the future. The proposed controller includes a few heuristic processes; rigorous certification of the control strategy is required, and it sometimes results in an extremely long controlling time to the target states. We believe that the concept of this controller explained in this study contributes to significantly decreasing number of thrusters in practical underactuated satellites.

The remainder of this paper is organized as follows. In Section 2, the formulations of the kinematic and dynamic equations of a FFRB are discussed. The procedure consists of a combination of the three hierarchical feedback steps and motion planning. Furthermore, the effects of the design parameters in each step are discussed. In Section 4, numerical simulations are conducted to verify the proposed controller, and a summary of the study is presented in Section 5.

2. Three-Dimensional Model of a FFRB (Free-Floating Rigid Body with Three Force Inputs)

2.1. Thruster Configuration. In order to simplify the discussion in Section 3, this subsection defines the three-thruster configuration for controlling a FFRB. It should be noted that the main purpose of this study is to demonstrate that a range of three-thruster configurations controls both the position and the attitude of a FFRB to the target values.

It is assumed that a rigid body has the three thrusters that satisfy the following conditions (i)-(iii) for its torques.

(i) Bidirectional control torques are generated about principal axes I and II of the FFRB.

This condition implies that the attitude motion in three dimensions is controlled by applying nonlinear control methods utilizing nonholonomic constraint of a FFRB unless the inertial moments of principal axes I and II are equal. (It is widely known that the attitude motion around the principal axis III is uncontrollable when the inertial moments are equal.)

(ii) With respect to the principal axis III, a unidirectional control torque is generated whenever control torques are applied around the principal axes I and II.

This condition simplifies the discussion in Section 3. Furthermore, with respect to the translational motion control, the following condition is assumed.

(iii) In the thruster configuration, the induced force along the principal axis III is bidirectional.

First, it is proven below that a range of three-thruster configurations satisfy the conditions (i) and (ii). We consider a thruster placed at [r.sub.i] and oriented to [d.sub.i] with respect to the body frame (see Figure 1(a)), and it generates a torque moment [T.sub.I] = [r.sub.i] x [d.sub.i] x [f.sub.i] with respect to the mass center. Three thrusters form a tetrahedron with their moments, and a range of three-thruster configurations have a projection of the tetrahedron on the I-II plane including the mass center inside the body as shown in Figure 1(b). Subsequently, the combination of the three forces in positive or zero values generates their resultant torque within the extension of the tetrahedron. Thus, this type of three-thruster configurations satisfies the conditions (i) and (ii). It is noted that a vertex of the tetrahedron is the mass center of the FFRB, and therefore, it is impossible to generate bidirectional control torques for all three principal axes with three thrusters as shown in [11, 12].

Next, it is shown that there is a three-thruster configuration that satisfies the condition (iii) under the conditions (i) and (ii) as follows. As described above, the moment vector [T.sub.i] is defined as the outer product between the thruster's position vector [r.sub.i] and orientation vector [d.sub.i]. However, as shown in Figure 1(c), there are numerous combinations for the two vectors in the plane perpendicular to the vector [T.sub.i]. This implies that the condition (iii) is satisfied as well as the conditions (i) and (ii) when one of [d.sub.i] (i = 1~3) has a negative component along the axis III (as shown in Figure 1(a)) and one of the others has a positive component (as shown in Figure 1(c)). Figure 1(d) illustrates an example of this configuration. It should be noted the three-thruster configuration used in the simulations in this study satisfies the three conditions.

2.2. Attitude Equations. We consider a FFRB with constant moments of inertia and with the absence of external disturbance forces. Euler equations of the FFRB are described with its principal axes as follows:

[mathematical expression not reproducible], (1)


[mathematical expression not reproducible]. (2)

The attitude is defined by the Euler parameter q as follows:

[mathematical expression not reproducible]. (3)

The differentiation of q is expressed with the angular velocity [omega] as follows:

[??] = Q(q) [omega], (4)


[mathematical expression not reproducible]. (5)

The total torque input T is given by the summation of [T.sub.i] as follows:

[mathematical expression not reproducible], (6)

where f = [[[f.sub.1] [f.sub.2] [f.sub.3]].sup.T], B = [[b.sub.1] [b.sub.2] [b.sub.3]], and [b.sub.i] = [r.sub.i] x [d.sub.i].

2.3. Translational Equations. The translational equation of the FFRB is expressed as follows:

M[??] = F. (7)

In (7), the translational force F induced by the thrusters is related to the control torque T. For this purpose, we first consider the force vector F' that is expressed with respect to the body fixed frame. Given that [F'.sub.i] = [d.sub.i] x [f.sub.i], the total translational force vector is expressed as follows:

F' = [3.summation over (i=1)] [F'.sub.i] = [3.summation over (i=1)] [d.sub.i] x [f.sub.i] = Df, (8)

where D = [[d.sub.1] [d.sub.2] [d.sub.3]]. The matrix B in (6) is nonsingular when the three control torques are independent, and thus, the force vector f is formulated as [B.sup.-1]T. Thus, the transformation matrix A(q) is used, and the translational force F in the inertial frame is expressed with the control torque T as follows:

F = A(q)D[B.sup.-1]T, (9)


[mathematical expression not reproducible]. (10)

Equation (9) indicates that the translational motion of the considering FFRB is induced dependently with the rotational motion.

3. Controller Design

3.1. Strategy. Without loss of generality, the target state can be set as q = [[000 [+ or -] 1].sup.T] and p = [[0 0 0].sup.T]. Furthermore, at the target attitude, the principal axis III of the FFRB coincides with the Z-axis of the inertial flame.

The switching controller proposed in this study uses the following three invariant manifolds: [S.sub.1] = 0, [S.sub.2] = 0, and [S.sub.3] = 0. The first one, namely [S.sub.1] = 0, corresponds to the manifold and indicates that the FFRB rotates about its principal axis III that coincides with the Z-axis. Hereafter, we refer to this state as "single spinning." The values of [S.sub.2] and [S.sub.3] indicate the remained velocity and position errors, respectively, along the Z-axis when the single spinning stops. Figure 2 shows the schematic image of the control process. It should be noted that maintaining a single spinning ([S.sub.1] = 0) under [T.sub.I] = 0 and [T.sub.II] = 0 after achieving [S.sub.2] = 0 and [S.sub.3] = 0 causes [S.sub.2] and [S.sub.3] to be invariant (this property is explained in the following section). Thus, various time histories of torque [T.sub.III] are available to stop the single spinning at the target attitude while maintaining [S.sub.2] = 0 and [S.sub.3] = 0. Proper planning of the [T.sub.III] history controls the position and velocity along X and Y to the target.

Figure 3 schematically illustrates the switching controller. The controller for achieving [S.sub.1] = 0, [S.sub.2] = 0, and [S.sub.3] = 0 is composed of the following three hierarchical control steps: Step A, Step B, and Step C. In Step A, [S.sub.1] = 0 (i.e., [q.sub.I], [q.sub.II], [[omega].sub.I], and [[omega].sub.II] = 0) is achieved using feedback law. In Step B, intentional disturbances and, in Step A, controls are repeated such that [S.sub.2] = 0 is satisfied. Thus, [S.sub.1] = 0 and [S.sub.2] = 0 are satisfied at the end of Step B. Step C achieves [S.sub.1] = 0, [S.sub.2] = 0, and [S.sub.3] = 0 by disturbing [S.sub.1] and [S.sub.2] to nonzero values and by repeating Step B to satisfy [S.sub.1] = 0 and [S.sub.2] = 0 again. When [S.sub.1] = 0, [S.sub.2] = 0, and [S.sub.3] = 0 are satisfied, Step D finally controls the single spinning angle to the target angle by the time profile of [T.sub.III], which is planned beforehand such that the position (X and Y) in the inertial frame coincides with the target positions. These steps are explained individually below given the assumption that the resultant torque component along the principal axis III is always positive, i.e., the three-thruster configuration satisfies [T.sub.III] > 0 as shown in Figure 1(b).

3.2. Step A: Control to the Single Spinning about the Axis III. Define [S.sub.1] as follows:

[S.sub.1] = [q.sup.2.sub.I] + [q.sup.2.sub.II] + [[omega].sup.2.sub.I] + [[omega].sup.2.sub.II].

This implies that at [S.sub.1] = 0, only [q.sup.III], [q.sub.IV], and [[omega].sub.III] correspond to nonzero values. Thus, achieving [S.sub.1] = 0 implies that the FFRB motion is led into a spinning motion around the axis III, and the axis coincides with the Z-axis of the inertial frame.

The following control law for [T.sub.I] and [T.sub.II] as proposed in [11] is adopted to achieve [S.sub.1] = 0 in the present study.

[T.sub.I] = -[K.sub.I][[omega].sub.I] - [J.sub.I][g.sub.I]/2 - [J.sub.I][bar.K][[omega].sub.II], (12)

[T.sub.I] = -[K.sub.II][[omega].sub.II] - [J.sub.II][g.sub.II]/2 - [J.sub.II][bar.K][[omega].sub.I], (13)


[g.sub.I] = [q.sub.I][q.sub.IV] + [q.sub.II][q.sub.III] + 2 [J.sub.II] - [J.sub.III]/[J.sub.I] [[omega].sub.II] [[omega].sub.II], (14)

[g.sub.II] = [q.sub.II][q.sub.IV] + [q.sub.I][q.sub.III] + 2 [J.sub.III] - [J.sub.I]/[J.sub.II] [[omega].sub.III] [[omega].sub.I]. (15)

[K.sub.I] and [K.sub.II] are arbitrary positive gains. [bar.K] is determined to provide a negative [[omega].sub.III] anytime such that positive [T.sub.III] reduces the single spinning rate. Reference [11] specifies the [bar.K] as follows:

[mathematical expression not reproducible], (16)


[mathematical expression not reproducible]. (17)

It should be noted that when [[omega].sup.2.sub.I]-[[omega].sup.2.sub.II] approaches zero, the magnitude of [bar.K] approaches infinity, and subsequently, the magnitudes of [T.sub.I] and [T.sub.II] become extremely high. Thus, in order to avoid this singularity, the control inputs in (12) and (13) should be designed to saturate when [[omega].sup.2.sub.I]-[[omega].sup.2.sub.II becomes quite low. Furthermore, by setting [bar.K] = 0, the control torques [T.sub.I] and [T.sub.II] are used as disturbances. The reason is explained in Step B.

3.3. Step B: Control for the Translational Velocity in the Z -Direction. Step B begins after the FFRB is in single spinning, and thus, the principal axis III coincides with the Z-axis in the inertial frame and the control torques around the axes I and II are [T.sub.I] = [T.sub.II] = 0. Therefore, the third component of (7) is expressed from (9) as follows:

M[??] = [[gamma].sub.3][T.sub.III], (18)

where [[[[gamma].sub.1] [[gamma].sub.2] [[gamma].sub.3]].sup.T] = [DB.sup.-1] [[0 0 1].sup.T] denotes the effect of the control torque [T.sub.III] on the translational motion in the single spinning. Furthermore, from (1) in the single spinning, (18) is reexpressed as follows:

[mathematical expression not reproducible]. (19)

The above equation implies that the translational motion in the Z-direction depends on the time profile of the angular acceleration around the axis III. When the single spinning ceases, (19) is analytically integrated to derive the velocity variation [DELTA][??] as follows:

[mathematical expression not reproducible]. (20)

Therefore, the translation velocity along the Z-axis in the inertial frame when the single spinning ceases is expressed as follows:

[S.sub.2] = [??] - [[gamma].sub.3]/M [J.sub.III][[omega].sub.III]. (21)

Conversely, we consider the time differential of [S.sub.2] when a control torque vector T = [[0, 0, [T.sub.III]].sup.T] is applied to the FFRB in the single spinning. It is expressed from the third components of (1), (7), and (9) as follows:

[mathematical expression not reproducible]. (22)

This implies that the variation of [S.sub.2] is generated by [T.sub.I], [T.sup.II], [[omega].sub.I], and [[omega].sub.I]. Thus, the value of [S.sub.2] does not change when the FFRB's single spinning is maintained. Conversely, a disturbance that makes [S.sub.1] [not equal to] 0 is necessary to change the value of [S.sub.2]. Thus, Step B requires an intentional disturbance that results in a nonzero value of [S.sub.1] and the process of Step A for returning [S.sub.1] to zero. The value of [S.sub.2] should approach zero during a series of these processes.

Intentional disturbances that satisfy the above process are shown as follows. Given that [T.sub.III] > 0 is assumed in this section, the angular rate [[omega].sub.III] around the axis III should approach zero and satisfy [[omega].sub.III] < 0. This implies that time integration of [[??].sub.III] must be negative and from the third component of (1), the following expression is obtained:

[[??].sub.III] = [T.sub.III]/[J.sub.III] + [[sigma].sub.III] [[omega].sub.I] [[omega].sub.II]. (23)

The second term ([[sigma].sub.III] [[omega].sub.I] [[omega].sub.II]) in the right-hand side of the above equation should be negative for a high portion of the trajectory for [S.sub.1] [not equal to] 0 and [S.sub.1] = 0. Therefore, when [J.sub.I] > [J.sub.II] is assumed for the FFRB, the trajectory shown in Figure 4 on the [[omega].sub.I]-[[omega].sub.II] plane is a solution for the intentional disturbance and mainly in the second or fourth quadrant. However, [[sigma].sub.III] [[omega].sub.I] [[omega].sub.II] is not analytically integrable, and thus, it is not possible to determine whether this trajectory always provides a negative value in the integration of [[sigma].sub.III] [[omega].sub.I] [[omega].sub.II]. (In the simulation of Section 4, the negative integration of [[sigma].sub.III] [[omega].sub.I] [[omega].sub.II] of the trajectory for the designed intentional disturbances is numerically confirmed.) It should be noted that the rotational direction (clockwise or counterclockwise) of the designed trajectory in the [[omega].sub.I]-[[omega].sub.II] plane dictates the sign of [DELTA][S.sub.2] where [DELTA][S.sub.2] denotes a variation of [S.sub.2] in a sequence [S.sub.1] = 0 [right arrow] [S.sub.1] [not equal to] 0 [right arrow] [S.sub.1] = 0. Thus, according to the value of (21), the direction to pull [S.sub.2] towards zero is specified. Additionally, lower disturbance trajectories result in lower variations for [S.sub.2]. Thus, by adjusting the direction and magnitude of intentional disturbance trajectories, [S.sub.2] converges to zero without overshooting.

To derive the above intentional disturbance, law as a general form requires rigorous mathematical analysis. Thus, as an alternative, we preliminarily determine the constants, [D.sub.Ib], [D.sub.IIb], and [D.sub.S2], in the following disturbance laws of (24), (25), and (26), respectively, by using numerical simulations as follows:

[T.sub.I] = [D.sub.Ib] sgn ([S.sub.2])/[absolute value of [[OMEGA].sub.III]], (24)

[T.sub.II] = [D.sub.IIb] sgn ([S.sub.2])/[absolute value of [[OMEGA].sub.III]], (25)

for the time duration corresponding to the following:

[t.sub.S2] = [D.sub.s2] [absolute value of [S.sub.2]], (26)

where [[OMEGA].sub.III] denotes the angular velocity around the axis III when intentional disturbances are applied. The disturbances are divided by the magnitude of [[OMEGA].sub.III] to decrease the effect on [DELTA][S.sub.2] (the variation in [S.sub.2] is proportional to the magnitude of [[OMEGA].sub.III] owing to the gyro effect). Equation (26) denotes that the intentional disturbance decreases with the magnitude of [S.sub.2] under positive [D.sub.S2]. The calculation of [DELTA][S.sub.2] requires numerical integration in [S.sub.1] [not equal to] 0, and thus, the setting process of [D.sub.Ib], [D.sub.IIb], and [D.sub.S2] is heuristic. Therefore, the saturation of [t.sub.S2] is desirable to avoid unexpected [DELTA][S.sub.2] behavior preliminary calculated by numerical simulations.

In order to cancel the effect of the [[OMEGA].sub.III] derivation, the following simple feedback law backs [[omega].sub.III] to the initial [[OMEGA].sub.III] as follows:

[T.sub.III] = -[K.sub.III] ([[omega].sub.III] - [[OMEGA].sub.III]), (27)

where [K.sub.III] denotes the positive gain. As shown in (22), (27) does not affect [S.sub.2].

3.4. Step C: Control of the Position in the Z-Direction. Step C controls the FFRB's position in the Z-direction in addition to Step B by achieving [S.sub.1] = 0, [S.sub.2] = 0, and [S.sub.3] = 0. In this step, [S.sub.3] is defined as the FFRB position when the single spinning ceases after the maneuver in Step D. As given below, [theta] = 2n[pi] is assumed as an identified preferable attitude angle in Step D where n denotes an arbitrary integer that indicates a rotation number in the single spinning.

The rotational angle [theta] about axis III is expressed as follows:

[theta] = 2 [tan.sup.-1] ([q.sub.III]/[q.sub.IV]) + 2n[pi]. (28)

If the stopping maneuver for the spinning begins at [theta] = [[theta].sub.3] + 2n[pi] and ends at [theta] = 2n[pi] in the single spinning, the variation of the Z-directional position, [DELTA]Z, is derived as follows:

[mathematical expression not reproducible], (29)

where [[tau].sub.d] denotes the time of drift motion before the final stopping maneuver in Step D. From [S.sub.2] = 0 in Step C, the drift interval is calculated from (21) as follows:

[[tau].sub.d] = [[theta].sub.3] + 2n[pi] - [theta]/[[omega].sub.III] = [J.sub.III][[gamma].sub.3]/M [[theta].sub.3] + 2n[pi] - [theta]/[??]. (30)

Substituting (30) into (29) yields the following expression for [DELTA]Z:

[DELTA]Z = [J.sub.III][[gamma].sub.3]/M (-[theta] + 2n[pi]). (31)

Thus, [S.sub.3] is defined as follows:

[S.sub.3] = Z + [DELTA]Z = [[xi].sub.1] + [[xi].sub.2]n, (32)


[[xi].sub.1] = Z - [J.sub.III][[gamma].sub.3]/M [theta], [[xi] = 2[pi][J.sub.III][[gamma].sub.3]/M. (33)

It should be noted that [S.sub.3] does not include [[theta].sub.3]. Subsequently, the time derivative of [S.sub.3] is derived as follows:

[mathematical expression not reproducible]. (34)

This equation indicates that [S.sub.3] is invariant under [S.sub.2] = 0 ([S.sub.1] = 0 is also satisfied by the control in Step B).

In (32), the value of [S.sub.3] depends on n. However, specifying a number for n in advance is not possible, because the position variation is uncertain until [S.sub.3] = 0 is achieved. In order to solve this problem, [S.sub.3] =0 are divided by [[xi].sub.2] to obtain [[xi].sub.1]/[[xi].sub.2] + n = 0. Subsequently, in order to ensure that n is an integer, [S.sub.3] is modified without n as follows:

[mathematical expression not reproducible], (35)

where MOD indicates the remainder. It should be noted that the modified [S.sub.3] is also invariant under [S.sub.2] = 0, because [[xi].sub.2] is a constant value. Furthermore, (35) has two forms to avoid discontinuity for the feedback control because the term MOD([[xi].sub.1], [[xi].sub.2]) is not continuous at MOD ([[xi].sub.1], [[xi].sub.2]) = 0. This discontinuity significantly affects the convergence of [S.sub.3] = 0.

Given that [S.sub.3] is invariant under [S.sub.2] = 0, the repetition of the sequence [S.sub.2] = 0 [right arrow] [S.sub.2] [not equal to] 0 [right arrow] [S.sub.2] = 0 is necessary to achieve [S.sub.3] = 0. The design procedure of the intentional disturbance for [S.sub.2] [not equal to] 0 is almost identical to that for [S.sub.1] [not equal to] 0 in Step B. The direction of the disturbance must be in the second or fourth quadrant on the [[[omega].sub.I]-[[omega].sub.II] plane in the case of [J.sub.I] > [J.sub.II], and the magnitude of the disturbance should be adjusted to avoid generating an excessive variation in [DELTA][S.sub.3] in a sequence.

Additionally, in Step B, we preliminary determine [D.sub.Ic], [D.sub.IIc], and [D.sub.S3] in the following disturbance law of (36), (37), and (38), respectively, by using numerical simulations as follows:

[T.sub.I] = [D.sub.Ic] sgn ([S.sub.3])/[absolute value of [[OMEGA].sub.III], (36)

[T.sub.II] = [D.sub.IIc] sgn ([S.sub.3]/[absolute value of [[OMEGA].sub.III], (37)

for the following time duration:

[t.sub.s3] = [D.sub.s3] [absolute value of [S.sub.3]]. (38)

3.5. Step D: Reaching the Target State. In Steps A, B, and C, [S.sub.1] =0, [S.sub.2] = 0, and [S.sub.3] = 0 were achieved, respectively. This means that the attitude angles around the axes I/II and the velocity/position along the Z-direction were adjusted for a target state. Thus, the remaining states correspond to the attitude angle around the axis III and the FFRB positions along X- and Y-directions. The remaining attitude angle is easily controlled by [T.sub.III], and thus, Step D adjusts the FFRB velocities/positions along the X- and Y-directions to the target positions by using a properly designed [T.sub.III] time profile until the single spinning ceases. Although there are various proper histories for [T.sub.III], this section proposes a solution for the motion planning problem.

The proposed motion planning procedure is composed of a series of control sequences. In each sequence, while [T.sub.I] = [T.sub.II] = 0 to maintain the single spinning state, the following [T.sub.III] input is used:

[T.sub.III] = [eta][??], (39)

where [eta] denotes a constant and is selected to satisfy [T.sub.III] > 0. It is assumed that the control maneuver is applied from t = [t.sub.0] to t = [t.sub.f]. In order to determine the proper [eta] and the time duration, the variations of the velocity and position along the X- and Y-directions are calculated as follows. First, from (7) and (9) for the control torque in (39) with [T.sub.I] = [T.sub.II] = 0, the velocity variation [mathematical expression not reproducible], generated in the j-th control maneuver is calculated as follows:

[mathematical expression not reproducible], (40)

where the subscripts 0 and f represent the values at the beginning and final instance of the j-th control maneuver, respectively. Subsequently, the position variation [mathematical expression not reproducible], is obtained from the integration of (40) as follows:

[mathematical expression not reproducible], (41)

where [DELTA][t.sub.j] denotes the time interval of this maneuver and [mathematical expression not reproducible] denotes the velocity vector at the beginning of the maneuver. With respect to the single spinning control with (39), the interval [DELTA][t.sub.j] is described as follows:

[DELTA][t.sub.j] = [J.sub.III]/[[eta].sub.j] ln (1 + [[theta].sub.jf] - [[theta].sub.j0]/[psi].sub.j]), (42)


[[phi].sub.j] = [J.sub.III][[??].sub.j0]/[[eta].sub.j]. (43)

Furthermore, in (41), [S.sub.[theta]]([[theta].sub.j0], [[theta].sub.jf], [[??].sub.j0], [[eta].sub.j]) and [C.sub.[theta]] ([[theta].sub.j0], [[theta].sub.jf], [[??].sub.j0], [[eta].sub.j]) are defined as the time integrations of sin ([theta]) and cos ([theta]) multiplied by [[eta].sub.J]/[J.sub.III], respectively, and expressed as follows:

[mathematical expression not reproducible], (44)

[mathematical expression not reproducible], (45)

where Si(*) and Ci(*) denote the sine and cosine integral functions, respectively, which are defined as follows:

[mathematical expression not reproducible]. (46)

The angular velocity [[??].sub.j+1)0] at the beginning of the (j + 1)-th maneuver in (43) is expressed in an algebraic form with [[eta].sub.j], [[theta].sub.j0], and [[theta].sub.jf] as follows:

[mathematical expression not reproducible], (47)

where [mathematical expression not reproducible], are determined by [[eta].sub.j], [[theta].sub.j0], and [[theta].sub.jf]. In this study, three control sequences (j = 1-3) are used to design the history of [theta]. The schematic image of the three sequences is shown in Figure 5. The details to specify the control parameters in each sequence are as follows:

The first sequence causes the velocities on the X-Y plane to correspond to zero. With respect to this control, [[eta].sub.1], [[theta].sub.10], and [[theta].sub.1f] are determined from (32) such that the following equation is satisfied:

[mathematical expression not reproducible], (48)

where [[theta].sub.10] and [[theta].sub.1f] are specified based on the direction of [[??].sub.10] and [[eta].sub.1 is determined from the norm of [[??].sub.10]. Following the maneuver, the FFRB is in rotational drift motion with [T.sub.III] = 0 until the second sequence begins.

Next, the combination of the second and the third sequences achieves [??] = 0. In the second sequence, with respect to [[theta].sub.2f] satisfying [[theta].sub.2f]-[[theta].sub.20] = -2[pi], (32) indicates that [DELTA][[??].sub.2] corresponds to zero (thus, [[??].sub.30] also corresponds to zero). Subsequently, given that [mathematical expression not reproducible] is expressed as follows:

[mathematical expression not reproducible], (49)


[mathematical expression not reproducible], (50)

[mathematical expression not reproducible], (51)

where [phi]([[theta].sub.j0], [[theta].sub.jf], [[??].sub.j0], [[eta].sub.j]) denotes the direction of the position variation. From (44), (45), (50), and (51), the value of [[eta].sub.2] affects the norm of [DELTA][[??].sub.2] and the value of [[theta].sub.20 affects the direction of [DELTA][[??].sub.2]. Prior to detailing the procedure to determine [[eta].sub.2] and [[theta].sub.20], we explain the third sequence because the value of [DELTA][[??].sub.3] is required to determine [[eta].sub.2] and [[theta].sub.20]. The third sequence begins when the drift motion brings the FFRB towards Z = 0 within a full rotation of [theta]. In the third sequence, the interval [[theta].sub.3f]-[[theta].sub.30] is set as -2[pi] for [DELTA][[??].sub.3] = 0 and [[theta].sub.3f] must correspond to 2(n + 1)[pi] to achieve the target attitude when the single spinning is stopped. Thus, [[eta].sub.3] is specified from (47) as follows:

[[eta].sub.3] = [J.sub.III][[??].sub.30]/2[pi], (52)

and [DELTA][[??].sub.3] is expressed as follows:

[mathematical expression not reproducible]. (53)

Thus, substituting (52) into (43) derives [[phi].sub.3] = -2[pi]. This indicates that [S.sub.[theta]] (2[pi], 0, [[??].sub.30], [[eta].sub.3]), [C.sub.[theta]] (2[pi], 0, [[??].sub.30], [[eta].sub.3]), and [phi](2[pi], 0, [[??].sub.30], [[eta].sub.3]) in (53) are not dependent on [[eta].sub.3] and [[??].sub.30], and thus, [DELTA][[??].sub.3] is a constant vector specified by a FFRB model and not by control maneuvers. Therefore, [[eta].sub.2] and [[theta].sub.20] are determined as follows. In order to obtain [??] = 0, [DELTA][[??].sub.2] + [DELTA][[??].sub.3] must correspond to -[[??].sub.20]. Subsequently, [DELTA][[??].sub.3] is a constant vector, and thus, [DELTA][[??].sub.2] is adjusted to satisfy [DELTA][[??].sub.2] + [DELTA][[??].sub.3] = -[[??].sub.20]. Consequently, [[eta].sub.2] is numerically calculated from (50), and (49) and (51) then yield the value of [DELTA][[??].sub.20] for the required [DELTA][[??].sub.2].

3.6. Final Form of the Controller with a Supplemental Sequence. It should be noted that [absolute value of [DELTA][[??].sub.2]] has an upper limit that depends on [[??].sub.20]. When a required [absolute value of [DELTA][[??].sub.2]] exceeds the upper limit, Step D does not achieve the target state. In order to avoid such a situation, minimizing [??] + [DELTA][[??].sub.3] is mandatory prior to Step D.

An idea for this minimization involves adjusting the timing when a sequence of the intentional disturbance and Step A ([S.sub.1] = 0 [right arrow] [S.sub.1] [not equal to] 0 [right arrow] [S.sub.1] = 0) begins. As shown above, the attitude angle 0 does not affect the controls of Steps A, B, and C. Thus, we switch on the control sequences while numerically predicting [mathematical expression not reproducible] to minimize the following V:

[mathematical expression not reproducible], (54)

where P and [??] denote the positions and velocities, respectively, immediately after each sequence [S.sub.1] = 0 [right arrow] [S.sub.1] [not equal to] 0 [right arrow] [S.sub.1] = 0 and V' denotes the squared distance between [??] = -[DELTA][[??].sub.3] and [??] = P. Numerical calculations are performed during the drift motions. The value L was heuristically determined given the following consideration. With respect to V' [greater than or equal to] L, V denotes the time derivative of [[absolute value of P + [DELTA][[??].sub.3]].sup.2] after each sequence. Thus, it implies that the maneuver begins at the attitude angle that minimizes the velocity vector for [??] = -[DELTA][[??].sub.3]. With respect to V' < L, the maneuver begins at the attitude angle minimizing the norm of [??]. For reference, we also attempted to minimize [mathematical expression not reproducible], which is a unified form of (54).

We emphasize that this supplement sequence is not always successful. Specifically, L requires severe tuning through trial and error in numerical simulations. This is the weak point in our proposed controller, and it is important to solve this problem in a future study.

The diagram of these sequences is shown as the "final control form" in Figure 6 as a summary of all the control sequences.

4. Example of the Proposed Controller with Numerical Simulation

This section describes the numerical validation of the control method shown in Section 3. The FFRB model and the controllers used in the simulations are summarized below, and the results of the numerical simulations are then shown.

4.1. FFRB Model and Initial State Value. According to the basic assumption of the study, the FFRB is assumed as with a constant mass, inertia, and no-upper-limitation thruster forces. In this example, the inertial moments of the FFRB are assumed as [J.sub.I] = 15 kg x [m.sup.2], [J.sub.II] = 10 kg x [m.sup.2], and [J.sup.III] = 20 kg x [m.sup.2]. Furthermore, the FFRB mass M is 20 kg. The placements and directions of the thrusters are set as follows:

Thruster 1: [r.sub.1] = [[0.0 1.0 - 1.0].sup.T], [d.sub.1] = [[0.0 0.0 - 1.0].sup.T]

Thruster 2: [r.sub.2] = [[-1.0 0.0 - 1.0].sup.T], [d.sub.2] = [[0.0 0.0 - 1.0].sup.T]

Thruster 3: [r.sub.3] = [[0.01 - 0.01 - 1.3].sup.T], [d.sub.3] = [[0.43 - 0.25 0.87].sup.T]

The forces of these thrusters are only positive, and their magnitudes are assumed as continuously changeable and without any upper limits. It should be noted again that the thruster configuration satisfies the conditions shown in Section 2.1. These model parameters result in [DELTA][[??].sub.3] = [[668 383].sup.T] for Step D.

Table 1 shows the initial state for the three cases used in simulations. Only the initial attitude is different between cases 1 and 2, and the initial position is different between cases 1 and 3. The angular and translational velocities are zero in all cases.

4.2. Controller

Step A. The parameters for the controller in (13) and (14) are set as [K.sub.1] =3.0, [K.sub.2] = 3.0, and [kappa] = 10, respectively.

Step B. The parameters for disturbance law in (24), (25), and (26) are set as [D.sub.Ib] = -0.3, [D.sub.IIb] = 0.3, and [D.sub.S2] = 0.5, respectively. Saturation is applied as [t.sub.S2] = 0.1 s for [absolute value of [S.sub.2]] > 0.1. Additionally, [K.sup.III] in (27) is set to 30.

Step C. The parameters for the disturbance law in (36), (37), and (38) are set as [D.sub.Ic] = 01, [D.sub.IIc] = -0.1, and [D.sub.S3] = 20, respectively. Saturation is applied as [t.sub.S3] = 0.1 s for [absolute value of [S.sub.3]] > 0.02.

Step D. As shown in the next subsection for simulations, the maneuver minimizing V that results from switch timing adjustments by using (54) leads to [??] = 0 after Step C in all the three cases. Thus, the first sequence to achieve zero velocities along the X- and Y-directions is omitted. The values for [[theta].sub.20] and [[eta].sub.2] specified are shown in the following subsection.

Supplemental sequence using (54): L in (54) is set to 200,000.

4.3. Simulation Results. Figures 7-10 show the results for case 1. As shown in Figures 7(a) and 7(b), the attitude successfully achieved the target state. Figure 7(a) indicates that the single spinning was dominant after initially achieving [S.sub.1] = 0. Figure 7(c) shows that the attitude velocities were also controlled to zero. Figure 8 verifies the success of the control of the positions. Figure 8 shows the time histories of [S.sub.1], [S.sub.2], and [S.sub.3] after initially achieving [S.sub.2] = 0. They confirm that the repetitions of the sequences, disturbances, and controls for [S.sub.1] = 0 and [S.sub.2] = 0 cause [S.sub.3] to quasi-asymptotically converge to zero. The thruster forces maintained positive values as shown in Figure 10, and the values of [[theta].sub.20] and [[eta].sub.2] were set as 3.27 rad and -0.174, respectively.

Cases 2 and 3 also show the success of the controls as illustrated in Figures 11-14. A significantly long period is involved in the process in case 2 owing to the low velocity in the Z-direction. This problem is caused by the small [[??]([approximately equal to] [[OMEGA].sub.III]), which lengthens the drifting time for approaching the target position, after the first Step A for [S.sub.1] = 0. The problem can be solved by adding acceleration step for [[OMEGA].sub.III] or tuning the parameters for Step A based on the initial states. The thruster forces maintained positive values in both cases 2 and 3 (the figures are omitted here). The values of [[theta].sub.20] and [[eta].sub.2] were set as 6.36 rad and -1.085, respectively, in case 2 and 4.55 rad and -0.179, respectively, in case 3.

5. Summary

This study presented a control method for three-dimensional attitude and position control of a free-floating rigid body using three thruster forces in which the directions are fixed to the body. First, the dynamics of the rigid body were described along with a thruster configuration to generate independent torques for two principal axes. The derived equations of motion for the rotation and translation confirmed that the equations were coupled. Next, a combination of switching controllers and motion planning was proposed. The switching controllers are composed of three hierarchical steps, and three invariant manifolds are used as switching surfaces. A process involving the application of an intentional disturbance and returning the state to on the lower class invariant manifolds is applied, and the states in upper class are individually controlled. In this control scheme, the magnitude of the intentional disturbance is tuned based on the distance to the invariant manifold. Motion planning in the final phase finalizes the state values to target values. The numerical simulation results verified the effectiveness of the proposed methods. In the simulations, the target attitudes and positions were successfully achieved in three cases with different initial states.

Although the limitations of this control method were not examined in detail, we first proposed a solution to this challenging problem. We believe that the methodology shown in this paper will contribute a practical satellite control in an underactuated situation due to thruster malfunction. A discussion of the generality of the method and overcoming the limitation will be addressed in the next phase of a future study.

a = [[a.sub.I] [a.sub.II]     Rotational axis with respect to
[a.sub.III].sup.T]:           the inertial frame

A(q) :                        Coordinate transformation
                              matrix from body coordinates
                              to the inertial frame

[bar.A] ([theta]) :           Rotation matrix in the X-Y

[d.sub.i]:                    Direction of the ith thruster
                              (i = 1-3) with respect to the
                              principal inertia axes of the
                              free-floating rigid body

[f.sub.i]:                    Magnitude of force generated by
                              the ith thruster (i = 1-3)

F:                            Total force vector generated by
                              the three thrusters with respect
                              to the inertial frame

F':                           Force vector generated by the
                              three thrusters with respect to
                              the body coordinates

[F'.sub.i]:                   Force vector generated by the
                              th thruster with respect to the
                              body coordinates

[J.sub.I],                    FFRB's principal moment of
[J.sub.II]                    inertia

M:                            Mass of the FFRB

m, n:                         Arbitrary integer

p = [[X Y Z].sup.T]           Position vector of the FFRB

[??] = [[X Y].sup.T]          Position vector of the FFRB on
                              the X-Y plane

[mathematical expression      Variations of [mathematical
not reproducible]:            expression not reproducible]
                              caused by the jth sequence in
                              Step D

q = [[q.sub.I]                Euler parameter for attitude
[q.sub.II] [q.sub.III]        definition with respect to the
[q.sub.IV]].sup.T]:           inertial frame

[r.sub.i]:                    Placement of the ith thruster
                              (i = 1-3) with respect to the
                              principal inertia axes of the

[DELTA][t.sub.j]:             Time duration for the maneuver
                              in the jth sequence of Step D

T = [[T.sub.I] [T.sub.II]     Total torque vector generated
[T.sub.III]].sup.T]:          by the thrusters around the
                              principal inertia axes of the

[T.sub.i]:                    Torque vector generated by the i
                              th thruster (i = 1-3) around the
                              principal inertia axes of the

[psi]:                        Rotational angle with respect to
                              [[[a.sub.I] [a.sub.II]

[theta]:                      Attitude angle with respect to
                              axis III during the single spin

[[tau].sub.d]:                Drift time before a stopping
                              maneuver for the single spin

[omega] = [[[omega].sub.I]    Angular velocity vector with
[[omega].sub.II]              respect to the principal inertia
[[omega].sub.III]].sup.T]:    axes of the FFRB

[[eta].sub.j],                Motion planning parameter in
[[phi].sub.j]                 the jth sequence of Step D

sgn ():                       Signum function

MOD (a, b):                   Remainder of alb

Subscription "j0" and "jf":   Initial and final values of the
                              jth sequence in Step D.

Data Availability

The authors confirm that the data supporting the findings of this study are available within the article.

Conflicts of Interest

The authors declare that they have no conflicts of interest.


[1] R. W. Brockett, R. S. Milman, and H. J. Sussmann, Eds., "Asymptotic stability and feedback stabilization," in Differential Geometric Control Theory, pp. 181-191, Birkhauser, Boston, USA, 1983.

[2] A. Isidori, Nonlinear Control Systems, Springer, 3rd edition, 1995.

[3] S. Hokamoto and M. Iwase, "Position and attitude control of a planar satellite by two thrusters," in AIAA Guidance, Navigation, and Control Conference, pp. 4206-4211, Chicago, USA, August 2009.

[4] P. Tsiotras and J. M. Longuski, "A new parameterization of the attitude kinematics," The Journal of the Astronautical Sciences, vol. 43, no. 3, pp. 243-262, 1995.

[5] P. Tsiotras and V. Doumtchenko, "Control of spacecraft subject to actuator failures state-of-the-art and open problems," Journal of Astronautical Sciences, vol. 48, no. 2, pp. 337-358, 2000.

[6] P. Morin and C. Samson, "Time-varying exponential stabilization of a rigid spacecraft with two control torques," IEEE Transactions on Automatic Control, vol. 42, no. 4, pp. 528534, 1997.

[7] S. Kasai, H. Kojima, and M. Satoh, "Spacecraft attitude maneuver using two single-gimbal control moment gyros," Acta Astronautica, vol. 84, pp. 88-98, 2013.

[8] A. Zavoli, G. De Matteis, F. Giulietti, and G. Avanzini, "Single-axis pointing of an underactuated spacecraft equipped with two reaction wheels," Journal of Guidance, Control, and Dynamics, vol. 40, no. 6, pp. 1465-1471, 2017.

[9] Y. Geng, D. Song, and R. Sun, "Inverse optimal stabilization of an underactuated spacecraft using two wheels," in 2016 8th International Conference on Intelligent Human-Machine Systems and Cybernetics (IHMSC), pp. 279-282, Hangzhou, China, August 2016.

[10] M. J. Sidi, Spacecraft Dynamics & Control: A Practical Engineering Approach, vol. 9, Cambridge University Press, 1997.

[11] T. Matsuno, Y. Yoshimura, and S. Hokamoto, "Geometric conditions of thrusters for 3D attitude control of a free-floating rigid spacecraft," Mathematics in Engineering, Science & Aerospace (MESA), vol. 5, no. 1, pp. 83-95, 2014.

[12] Y. Yoshimura, T. Matsuno, and S. Hokamoto, "Three dimensional attitude control of an underactuated satellite with thrusters," International Journal of Automation Technology, vol. 5, no. 6, pp. 892-899, 2011.

[13] H. Khennouf and C. Canudas De Wit, "On the construction of stabilizing discontinuous controllers for nonholonomic systems," IFAC Proceedings Volumes, vol. 28, no. 14, pp. 667-672, 1995.

[14] R. Mukherjee and M. Kamon, "Almost smooth time-invariant control of planar space multibody systems," IEEE Transactions on Robotics and Automation, vol. 15, no. 2, pp. 268-280, 1999.

[15] T. Matsuno and S. Hokamoto, "Position and attitude feedback controller of a planar satellite with two thrusters," in ICROSSICE International Joint Conference, pp. 4206-4211, Fukuoka, Japan, August 2009.

[16] T. Matsuno and S. Hokamoto, "Position and attitude control of a planar satellite by two constant force thrusters," in AIAA/AAS Astrodynamics Specialist Conference, Toronto, Ontario, Canada, August 2010.

[17] Y. Yoshimura, T. Matsuno, and S. Hokamoto, "Global trajectory design for position and attitude control of an underactuated satellite," Transactions of the Japan Society for Aeronautical and Space Sciences, vol. 59, no. 3, pp. 107-114, 2016.

[18] H. J. Sussmann, "A general theorem on local controllability," SIAM Journal on Control and Optimization, vol. 25, no. 1, pp. 158-194, 1987.

[19] B. Goodwine and J. Burdick, "Controllability with unilateral control inputs," in Proceedings of 35th IEEE Conference on Decision and Control, pp. 3394-3399, Kobe, Japan, December 1996.

Takashi Matsuno (iD), (1) Yasuhiro Yoshimura, (2) and Shinji Hokamoto (2)

(1) Department of Mechanical and Aerospace Engineering, Tottori, Tottori University, 4-101 Koyama-cho, Tottri 680-8552, Japan

(2) Department of Aeronautics and Astronautics, Fukuoka, Kyushu University, 744 Motooka, Nishi-ku, Fukuoka 819-0395, Japan

Correspondence should be addressed to Takashi Matsuno;

Received 22 March 2018; Accepted 28 August 2018; Published 29 October 2018

Academic Editor: Paolo Gasbarri

Caption: Figure 1: Three-thruster configuration for attitude and position control. (a) Direction of the torque vector. (b) Geometry formed with the three torque vector. (c) Another thruster placement. (d) Three-thruster configuration.

Caption: Figure 2: Schematic image of the [S.sub.1], [S.sub.2], and [S.sub.3].

Caption: Figure 3: Schematic of the control sequence.

Caption: Figure 4: Trajectory of [[omega].sub.I] and [[omega].sub.II] in Step B in the case of [J.sub.I] > [J.sub.II].

Caption: Figure 5: Schematic images of motion planning with the three sequences.

Caption: Figure 6: Diagram of the final control form.

Caption: Figure 7: Attitude in case 1. (a) Attitude trajectory to the target. (b) Attitude history around the target. (c) Time history of attitude velocity.

Caption: Figure 8: Position in case 1. (a) Position trajectory to the target. (b) Velocity histories.

Caption: Figure 9: Time history of the three invariant manifolds after first achieving [S.sub.2] = 0 in case 1.

Caption: Figure 10: Time history of thruster forces in case 1. Thruster 1 generates [f.sub.1], thruster 2 generates [f.sub.2], and thruster 3 generates [f.sub.3].

Caption: Figure 11: Attitude in case 2. (a) Attitude trajectory to the target. (b) Attitude history around the target.

Caption: Figure 12: Position in case 2. (a) Position trajectory to the target. (b) Velocity histories around the target.

Caption: Figure 13: Attitude in case 3. (a) Attitude trajectory to the target. (b) Attitude history around the target.

Caption: Figure 14: Position in case 3. (a) Position trajectory to the target. (b) Velocity histories around the target.
Table 1: Initial state values for the numerical simulations.

           [q.sub.I]   [q.sub.II]   [q.sub.III]   [q.sub.IV]

Case 1        0.4        -0.25          0.1          0.876
Case 2        0.8        -0.25          0.1          0.536
Case 3        0.4        -0.25          0.1          0.876

             X (m)       Y (m)       Z (m)

Case 1       -200         700        -350
Case 2       -200         700        -350
Case 3        0.0         0.0        -500
COPYRIGHT 2018 Hindawi Limited
No portion of this article can be reproduced without the express written permission from the copyright holder.
Copyright 2018 Gale, Cengage Learning. All rights reserved.

Article Details
Printer friendly Cite/link Email Feedback
Title Annotation:Research Article
Author:Matsuno, Takashi; Yoshimura, Yasuhiro; Hokamoto, Shinji
Publication:International Journal of Aerospace Engineering
Article Type:Report
Geographic Code:1USA
Date:Jan 1, 2018
Previous Article:Development of Design Methodology for a Small Solar-Powered Unmanned Aerial Vehicle.
Next Article:Characteristics of Orbit Determination with Short-Arc Observation by an Optical Tracking Network, OWL-Net.

Terms of use | Privacy policy | Copyright © 2020 Farlex, Inc. | Feedback | For webmasters