Journal of Lie Theory Volume 7 (1997) 269–278 C ◯1997 Heldermann Verlag

Stability and control in spacecraft dynamics

Mircea Puta

Communicated by W. A. F. Ruppert

Abstract. An optimal control problem for the spacecraft dynamics is discussed and some of its properties are pointed out.

1. Introduction

Recent work in nonlinear control has drawn attention to drift–free systems with fewer controls than state variables. These arise in problems of motion planning for wheeled robots subject to nonholonomic controls [12], [13], models of kinematic drift effects in space systems subject to appendage vibrations or articulations [5], [6], models of self-propulsion of paramecia at low Reynolds numbers [16] and autonomous underwater vehicle dynamics [8].

The goal of our paper is to discuss a similar problem for the spacecraft dynamics. The case with drift can be found in [11].

2. The spacecraft as a drift-free left invariant system on SO(3)

Let us consider a spacecraft free to move in ℝ^{3} . Let (b_{1}, b_{2}, b_{3}) be an orthonormal frame fixed on the body and let (r_{1}, r_{2}, r_{3}) define an inertial frame with the origin coincident with the origin of the body-fixed frame. Then we define X(t) ∈ SO(3) ( = the special orthogonal group, i.e., SO(3) = {A ∈ M_{3×3}(ℝ) | A^{t} · A = I_{3}, det(A) = 1}) such that r_{i} = X(t)b_{i}, i.e., X(t) determines the attitude of the spacecraft at time t. Let e_{1} = (1, 0, 0)^{T} , e_{2} = (0, 1, 0)^{T} and e_{3} = (0, 0, 1)^{T} . Define

ˆ : x = [ \begin{matrix}x_{1} \\ x_{2} \\ x_{3} \end{matrix} ] ∈ ℝ^{3} ↦ x̂ = [ \begin{matrix}0 & −x_{3} & x_{2} \\ x_{3} & 0 & −x_{1} \\ −x_{2} & x_{1} & 0 \end{matrix} ] ∈ 𝔰𝔬(3)

and

A_{i} = ê_{i}; i = 1, 2, 3.

ISSN 0949–5932 / $2.50 ^{_{C} ◯} Heldermann Verlag Then {A_{1}, A_{2}, A_{3}} is the standard basis for 𝔰𝔬(3) and X(t) satisfies:

Ẋ = XΩ̂; Ω̂ = ∑_{i=1}^{3} Ω_{i}(t)A_{i},  (2.1)

where Ω = (Ω_{1}, Ω_{2}, Ω_{3})^{T} is the angular velocity of the spacecraft in the bodyfixed coordinates. If we let

u_{i} = Ω_{i}; i = 1, 2, 3,

i.e., if we interpret the components of the angular velocity as our controls, then (2.1) takes the form

X = X (∑_{i=1}^{3} u_{i}A_{i}) .  (2.2)

We shall be most interested in the case when only two components of the angular velocity can be controlled. For example, if we can control the angular velocity about the b_{1} and b_{2} axes, then X(t) satisfies

Ẋ = X(u_{1}A_{1} + u_{2}A_{2}).

This realization of the spacecraft dynamics is due to Leonard [7].

Theorem 2.1. The system (2.3) is controllable and it is a single bracket one.

Let us observe that the above theorem tells us in fact that we can reorient the spacecraft as desired by controlling only two of the three angular velocity components (e.g. roll and pitch velocities).

Remark 2.1. Similar results hold for the controls about:

(i) b_{1} and b_{2} axes; (ii) b_{2} and b_{3} axes.

3. An optimal control problem for the spacecraft dynamics

Let J be the cost function given by

J(u_{1}, u_{2}) = \frac{1}{2} ∫_{0}^{t_{f}}[c_{1}u_{1}^{2}(t) + c_{2}u_{2}^{2}(t)]dt; c_{1} > 0, c_{2} > 0.  (3.1)

Then we can prove:

Theorem 3.1. The controls that minimize J and steer the system (2.3) from X = 0 at t = 0 to X = X_{f} at t = t_{f} are given by:

u_{1} = \frac{1}{c_{1}} P_{1}; u_{2} = \frac{1}{c_{2}} P_{2},

where the functions P_{i} are solutions of

\{ \begin{matrix}Ṗ_{1} = −\frac{1}{c_{2}} P_{2}P_{3} \\ Ṗ_{2} = \frac{1}{c_{1}} P_{1}P_{3} \\ Ṗ_{3} = (\frac{1}{c_{2}} − \frac{1}{c_{1}} ) P_{1}P_{2} . \end{matrix}  (3.2)

Proof. Simply apply Krishnaprasad’s theorem [4]. It follows that the optimal Hamiltonian is given by

h = \frac{1}{2c_{1}} P_{1}^{2} + \frac{1}{2c_{2}} P_{2}^{2}.  (3.3)

It is in fact the controlled Hamiltonian H given by

H = P_{1}u_{1} + P_{2}u_{2} − \frac{1}{2}(c_{1}u_{1}^{2} + c_{2}u_{2}^{2}),

which is reduced to (𝔰𝔬(3))^{⋆} via the Poisson reduction. Here (𝔰𝔬(3))^{⋆} is (𝔰𝔬(3))^{⋆} together with the minus Lie–Poisson structure

Π_{RB} = [ \begin{matrix}0 & −P_{3} & P_{2} \\ P_{3} & 0 & −P_{1} \\ −P_{2} & P_{1} & 0 \end{matrix} ] .  (3.4)

Then the optimal controls are given by

u_{1} = \frac{1}{c_{1}} P_{1}; u_{2} = \frac{1}{c_{2}} P_{2},

where the functions P_{i} are solutions of the reduced Hamilton’s equations (or momentum equations) given by

[Ṗ_{1}, Ṗ_{2}, Ṗ_{3}] = Π_{RB} · ∇h,

which are nothing else but the required equations.

Remark 3.1. The function C given by

C = \frac{1}{2}(P_{1}^{2} + P_{2}^{2} + P_{3}^{2})  (3.5)

is a Casimir of our configuration ((𝔰𝔬3))^{⋆}, Π_{RB}) ≃ (ℝ^{3}, Π_{RB}), i.e.,

(∇C)^{T} · Π_{RB} = 0.

Remark 3.2. The phase space curves of our system (3.2) are the intersections of the elliptic cylinders

\frac{P_{1}^{2}}{c_{1}} + \frac{P_{2}^{2}}{c_{2}} = 2h

with the spheres

P_{1}^{2} + P_{2}^{2} + P_{3}^{2} = 2C.

Theorem 3.2. The dynamics (3.2) is equivalent to the pendulum dynamics. Proof. Indeed, h is a constant of motion, so

\frac{P_{1}^{2}}{c_{1}} + \frac{P_{2}^{2}}{c_{2}} = l^{2}.

Let us take now

\{ \begin{matrix}P_{1} = l\sqrt{c_{1}} cos θ \\ P_{2} = l\sqrt{c_{2}} sin θ . \end{matrix}

Then

Ṗ_{1} = −\sqrt{c_{1}} sin θ · θ̇

= −l\sqrt{\frac{c_{1}}{c_{2}}} \sqrt{c_{2}} sin θ · θ̇

= −\sqrt{\frac{c_{1}}{c_{2}}} P_{2}θ̇,

or equivalently,

θ̇ = −\sqrt{\frac{c_{1}}{c_{2}}} \frac{Ṗ_{1}}{P_{2}}

= −\sqrt{\frac{c_{1}}{c_{2}}} (−\frac{P_{2}P_{3}}{P_{2}} ) \frac{1}{c_{2}}

= \frac{1}{\sqrt{c_{1}c_{2}}} P_{3}.

Differentiating again, we get

θ̈ = \frac{l^{2}}{\sqrt{c_{1}c_{2}}} (\frac{1}{c_{2}} − \frac{1}{c_{1}} ) c_{1}c_{2} sin 2θ,

or equivalently

θ̈ = \frac{l^{2}}{2} \sqrt{c_{1}c_{2}}(c_{1} − c_{2}) sin 2θ.  (3.6)

Thus, pendulum dynamics as required.

In the particular case c_{1} = c_{2} = 1 we refined Baillieul’s theorem [2], namely:

Theorem (Baillieul). The controls which minimize

J(u_{1}, u_{2}) = \frac{1}{2} ∫_{0}^{t_{f}}(u_{1}^{2} + u_{2}^{2})dt

and steer the system (2.3) from X = 0 at t = 0 to X = X_{f} at t = t_{f} are given by sinusoids. Proof. Indeed, let us take in (3.6) c_{1} = c_{2} = 1. Then θ = k_{1}t + k_{2} and the optimal controls are given via theorems 3.1 and 3.2 by

\begin{matrix}u_{1} = l cos(k_{1}t + k_{2}) \\ u_{2} = l sin(k_{1}t + k_{2}) \end{matrix}

as required.

Theorem 3.4. The system (3.2) may be realized as a Hamilton–Poisson system in an infinite number of different ways, i.e., there exist infinitely many different (in general nonisomorphic) Poisson structures on ℝ^{3} such that the system (3.2) is induced by an appropriate Hamiltonian . Proof. Indeed, to begin with, let us observe that our system can be put in an equivalent form:

Ṗ = ∇C × ∇h,

where P = [P_{1}, P_{2}, P_{3}]^{T} and C, h are respectively given by (3.5) and (3.3). Now, an easy computation shows us that the system (3.2) may be realized as a Hamilton–Poisson system with the phase space ℝ^{3} , the Poisson bracket {·, ·}_{ab} given by

{f, g}_{ab} = −∇C′ · (∇f × ∇g),

where a, b ∈ ℝ and

C′ = \frac{a}{2}C + \frac{b}{2}h,

and the Hamiltonian h′ defined by

h′ = cC + dh,

where c, d ∈ ℝ, ad − bc = 1.

Finally, we shall discuss the integrability of the equations (3.2) via elliptic functions. More precisely, we have

Theorem 3.5. The equations (3.2) may be explicitely integrated by elliptic functions. Proof. It is known that

P_{1}^{2}c_{2} + P_{2}^{2}c_{1} = l, l = 2hc_{1}c_{2}

and

P_{1}^{2} + P_{2}^{2} + P_{3}^{2} = 2C

are constants of motion. Then an easy computation shows us that

P_{2}^{2} = \frac{c_{2}}{c_{2} − c_{1}} [\frac{2Cc_{2} − l}{c_{2}} − P_{3}^{2}]

and

P_{1}^{2} = \frac{c_{1}}{c_{1} − c_{2}} [\frac{2Cc_{1} − l}{c_{1}} − P_{3}^{2}] .

Using now the third equation from (3.2) we get

(Ṗ_{3})^{2} = \frac{1}{c_{1}c_{2}} (P_{3}^{2} − \frac{2Cc_{2} − l}{c_{2}} ) (\frac{2Cc_{1} − l}{c_{1}} − P_{3}^{2}) ,

that is

t = ∫_{P_{3}(0)}^{P_{2}} \frac{dt}{\sqrt{\frac{1}{c_{1}c_{2}} (P_{3}^{2} − \frac{2Cc_{2}−l}{c_{2}}) (\frac{2Cc_{1}−l}{c_{1}} − P_{3}^{2})}} ,

which shows that P_{3} , and hence P_{1}, P_{2} are elliptic functions of time.

4. Numerical integration of the equation (3.2)

In this section we shall discuss the numerical integration of the equations (3.2) via the Lie–Trotter formula and the midpoint rule, and we shall point out some of their geometric properties.

To begin with, let us observe that the Hamiltonian vector field X_{h} splits as follows:

X_{h} = X_{h_{1}} + X_{h_{2}},

where

h_{1} = \frac{1}{2c_{1}} P_{1}^{2} h_{2} = \frac{1}{2c_{2}} P_{2}^{2}.

The integral curves of X_{h_{1}} and X_{h_{2}} are given by

P (t) = exp(tX_{h_{1}}) · P (0) = Φ_{1}(t, P (0))

and respectively,

P (t) = exp(tX_{h_{2}}) · P (0) = Φ_{2}(t, P (0)).

Now, following [10] (see also [15]), the Lie–Trotter formula gives rise to an explicit integrator of the equation (3.2), namely:

P^{k+1} = Φ_{1}(t, Φ_{2}(t, P^{k}))

or explicitely:

P_{1}^{k+1} = P_{1}^{k} cos \frac{P_{2}(0)}{c_{2}}t − P_{3}^{k} sin \frac{P_{2}(0)}{c_{2}}t

P_{2}^{k+1} = P_{1}^{k} sin \frac{P_{1}(0)}{c_{1}}t sin \frac{P_{2}(0)}{c_{2}}t + P_{2}^{k} cos \frac{P_{1}(0)}{c_{1}}t + P_{3}^{k} sin \frac{P_{1}(0)}{c_{1}}t cos \frac{P_{1}(0)}{c_{1}}t

P_{3}^{k+1} = P_{1}^{k} cos \frac{P_{1}(0)}{c_{1}}t sin \frac{P_{2}(0)}{c_{2}}t − P_{2}^{k} sin \frac{P_{1}(0)}{c_{1}}t + P_{3}^{k} cos \frac{P_{1}(0)}{c_{1}}t cos \frac{P_{2}(0)}{c_{2}}t

(4.1)

Some of its properties are sketched in the following theorem:

Theorem 4.1.

(i) The numerical integrator (4.1) preserves the Poisson structure (3.4). (ii) The numerical integrator (4.1) preserves the Casimirs of our config-

uration (ℝ^{3}, Π_{RB}).

(iii) Its restriction to each coadjoint orbit (C = \frac{1}{2}k_{2}, ω_{k} = \frac{1}{k} (P_{2}dP_{1} ∧ dP_{3} − P_{3}dP_{1} ∧ dP_{2} − P_{1}dP_{2} ∧ dP_{3})) gives rise to a symplectic integrator.

(iv) The numerical integrator (4.1) does not preserve the Hamiltonian h given by (3.3). Proof. The items (i) – (iii) hold because Φ_{1} and Φ_{2} are flows of some Hamiltonian vector fields, hence they are Poisson maps.

Item (iv) is essentially due to the fact that {h_{1}, h_{2}} ≠ 0.

An alternative way for the numerical integration of the equations (3.2) is given by the midpoint rule of Gauss–Legendre integrator. In our case it is given by the following implicit formula:

\{ \begin{matrix}P_{1}^{k+1} = P_{1}^{k} − \frac{h}{4c_{2}} (P_{2}^{k+1} + P_{2}^{k})(P_{3}^{k+1} + P_{3}^{k}) \\ P_{2}^{k+1} = P_{2}^{k} + \frac{h}{4c_{1}} (P_{1}^{k+1} + P_{1}^{k})(P_{3}^{k+1} + P_{3}^{k}) \\ P_{3}^{k+1} = P_{3}^{k} + \frac{h(c_{1} − c_{2})}{4c_{1}c_{2}}(P_{1}^{k+1} + P_{1}^{k})(P_{2}^{k+1} + P_{2}^{k}) \end{matrix}  (4.2)

Using now the same arguments as in [1] with obvious modifications we can prove

Theorem 4.2.

(i) The numerical integrator (4.3) preserves the Hamiltonian (3.3) and the Casimir (3.5).

(ii) The numerical integrator (4.3) does not preserve the Poisson structure (3.4).

(iii) The numerical integrator (4.3) does not preserve the symplectic structure:

ω_{k} = \frac{1}{k} [P_{2}dP_{1} ∧ dP_{3} − P_{3}dP_{1} ∧ dP_{2} − P_{1}dP_{2} ∧ dP_{3}],

hence its restriction to each coadjoint orbit (C = \frac{1}{2}k^{2}, ω_{k}) is not a symplectic integrator.

5. Stability

It is not hard to see that the equilibrium states of our system (3.2) are

e_{1} = (M, 0, 0); e_{2} = (0, M, 0); e_{3} = (0, 0, M),

where M ∈ ℝ. Now we shall discuss their nonlinear stability. Recall that an equilibrium point p is nonlinearly stable if trajectories starting close to p stay close to p. In other words, a neighbourhood of p must be flow invariant.

First consider the system linearized about c_{1} . Its eigenvalues are given by solutions of the equation

λ (λ^{2} − \frac{c_{1} − c_{2}}{c_{1}^{2}c_{2}} M^{2}) = 0.

If c_{1} > c_{2} , then a root of the characteristic polynomial has positive real part, thus e_{1} is unstable. If c_{1} < c_{2} , then the characteristic polynomial has two imaginary eigenvalues and one zero eigenvalue. Is the system stable? We shall prove that it is, via the energy–Casimir method.

Recall that the energy–Casimir method (see [3], [9] or [14]) requires finding a constant of motion for the system, say H, usually the energy, and a family of constants of motion C such that for some C ∈ C, C + H has a critical point at the equilibrium of interest. C ’s are often taken to be Casimirs. Definiteness of δ^{2}(H + C), the second variation of H + C and the critical point is sufficient to prove the stability if the phase space of the system is finite dimensional. We use this method to prove the following

Theorem 5.1. The equilibrium state e_{1} = (M, 0, 0) is nonlinear stable if c_{1} < c_{2} .

Proof. Consider the energy–Casimir function

h_{φ} = \frac{1}{2} (\frac{P_{1}^{2}}{c_{1}} + \frac{P_{2}^{2}}{c_{2}} ) + φ (\frac{1}{2}(P_{1}^{2} + P_{2}^{2} + P_{3}^{2})) ,

where φ : ℝ → ℝ is an arbitrary smooth real valued function defined on ℝ. Let φ′, φ′′ denote its first and second derivatives. Now, the first variation of h_{φ} is given by

δh_{φ} = \frac{P_{1}}{c_{1}} δP_{1} + \frac{P_{2}}{c_{2}} δP_{2} + φ′(·)(P_{1}δP_{1} + P_{2}δP_{2} + P_{3}δP_{3}).

This equals zero at the equilibrium of interest if and only if

φ′ (\frac{1}{2}M^{2}) = −\frac{1}{c_{1}} .  (5.1)

The second variation of h_{φ} at the equilibrium of interest is given via (5.1) by

δ^{2}h_{φ}(c_{1}) = φ′′ (\frac{1}{2}M^{2}) M^{2}(δP_{1})^{2} + \frac{c_{1} − c_{2}}{c_{1}c_{2}} (δP_{2})^{2} − \frac{1}{c_{1}} (δP_{3})^{2}.

Since c_{1} < c_{2}; c_{1}, c_{2} > 0 and having chosen φ such that

φ′′ (\frac{1}{2}M^{2}) < 0,

we can conclude that the second variation at the equilibrium of interest is negative definite and thus e_{1} is nonlinear stable.

In a similar manner we can prove

Theorem 5.2. The equilibrium state e_{2} = (0, M, 0) is (i) unstable, if c_{1} < c_{2} ;

(ii) nonlinear stable, if c_{1} > c_{2} .

Theorem 5.3. The equilibrium state e_{3} = (0, 0, M) is nonlinear stable.

References

[1] Austin, M. A., P. S. Krishnaprasad, and L. S. Wang, Almost Poisson integration of rigid body systems, J. Comp. Phys. 107 (1993), 105–117.

[2] Baillieul, J., “Some optimization problems in control theory,” Ph.D. Thesis, Harvard University, 1975.

[3] Holm, D., J. Marsden, T. Raţiu, and A. Weinstein, Nonlinear stability of fluid and plasma equilibria, Phys. Reports 123 (1985), 1–116.

[4] Krishnaprasad, P. S., “Optimal control and Poisson reduction,” Technical Report 93–87, Institute for System Research, University of Maryland, 1993.

[5] Krishnaprasad, P. S., and R. Yang, Geometric phases, anholonomy and optimal movement, in Proc. IEEE Int. Conf. Robot Automat. Sacramento CA (1991), 2185–2189.

[6] Krishnaprasad, P. S., R. Yang, and W. P. Dayawansa, Control problems on principal bundles and nonholonomic mechanics, in Proc. 30th IEEE Conf. Decis. Control. Brighton UK (1991), 1133–1138.

[7] Leonard, N. E., “Averaging and motion control of systems on Lie groups,” Ph.D. Thesis, University of Maryland, College Park, MD, 1994.

[8] Leonard N. E., and P. S. Krishnaprasad, Motion control on drift–free, left invariant systems on Lie groups, IEEE Trans. on Automatic Control 40 (1995), 1539–1554.

[9] Marsden, J., and T. Raţiu “Introduction to Mechanics and Symmetry,” Springer-Verlag, Berlin etc., 1994.

[10] McLachlan, R., On the numerical integration of ordinary differential equations by symmetric composition methods, SIAM J. Sci. Comp. 16 (1995), 151–168.

[11] Montgomery, R. S., S. Sastry, and G. Walsh, Optimal path planning on Lie groups, Proc. 1994 IEEE Conf. Decis. Contr. (1995), 1258–1263.

[12] Murray, R. M., and S. S. Sastry, Steering nonholonomic systems using sinusoids, in: Proc. 29th IEEE Conf. Decis. Contr., Honolulu (1990), 2097–2101.

[13] —, Nonholonomic motion planning: Steering using sinusoids, IEEE Trans. Automat. Contr. 38 (1993), 700–716.

[14] Puta, M. “Hamiltonian mechanical systems and geometric quantization,” Math. and Appl., 260, Kluwer Academic Publishers, 1993.

[15] —, Lie–Trotter formula and Poisson dynamics, J. Diff. Geom. in Balkan, to appear.

[16] Shapere, A., and F. Wilczek, Geometry of self-propulsion at low Reynolds number, J. Fluid Mech. 198 (1989), 557–585.

278 Puta

Seminarul de Geometrie–Topologie, University of Timişoara, B-dul V. Pârvan 4, 1900, Timişoara, Romania Fax: 0040-56-194002 puta@tim1.uvt.ro

Received August 8, 1996 and in final form April 14, 1997