Close
About
FAQ
Home
Collections
Login
USC Login
Register
0
Selected
Invert selection
Deselect all
Deselect all
Click here to refresh results
Click here to refresh results
USC
/
Digital Library
/
University of Southern California Dissertations and Theses
/
New approaches to satellite formation-keeping and the inverse problem of the calculus of variations
(USC Thesis Other)
New approaches to satellite formation-keeping and the inverse problem of the calculus of variations
PDF
Download
Share
Open document
Flip pages
Contact Us
Contact Us
Copy asset link
Request this asset
Transcript (if available)
Content
NEW APPROACHES TO SATELLITE FORMATION-KEEPING AND
THE INVERSE PROBLEM OF THE CALCULUS OF VARIATIONS
by
Hancheol Cho
A Dissertation Presented to the
FACULTY OF THE USC GRADUATE SCHOOL
UNIVERSITY OF SOUTHERN CALIFORNIA
In Partial Fulfillment of the
Requirements for the Degree
DOCTOR OF PHILOSOPHY
(AEROSPACE ENGINEERING)
December 2012
Copyright 2012 Hancheol Cho
ii
Acknowledgements
I would like to convey my special thank to my advisor Professor Firdaus E. Udwadia.
Without him, this dissertation may not be possible to be published. While I studied at the
University of Southern California, he was always willing to help, guide, support me so
that I could survive as an independent researcher. I especially learned from him the most
important things in this world: integrity, humility, courage, and helpfulness. I will always
cherish his words as long as I live; “We can change this world by helping other people,
not by delivering our knowledge.”
I am also deeply grateful to my dissertation committee members Professor Bingen Yang
and Professor Daniel Erwin who have made invaluable comments on my research and
provided a great deal of encouragement over the years. I am also deeply indebted to
Professor Henryk Flashner, Professor Paul Newton, Professor Michael Gruntman,
Professor Chunming Wang, and Dr. Gerald Hintz for their inspiration on my research
topics.
I also would like to greatly appreciate the people at the Power of Praise Church and Pastor
Synn. They always encourage me and pray for me when I have hard time in research or
human relations. I cannot recollect my USC years without them.
iii
Finally I want to say I really really love my family. My dad, mom, and younger brothers
(Hanmin and Hanjin) in South Korea always have supported my study in both financially
and spiritually. I am so sorry for having been careless about my family due to study, and
want to be a good son after completing my study.
I never forget the fact that God always leads my every step. I am still curious about what
God will display in my life. Now I ask God’s blessing on the great adventure I am going to
set sail.
iv
Table of Contents
Acknowledgements……………………………………………………….………..
List of Figures………………………………………………………………..……..
List of Tables……………………………………………………………………….
Abstract……………………………………………………………….….…………
Chapter 1 GENERAL INTRODUCTION…………………………………………
1.1 Motivation…………………………………………………...…………..
1.2 Background……………………………………………………….……...
1.3 Research Objectives……………………………………………….…….
1.4 The Study’s Contributions………………………………………….……
1.5 Organization……………………………………………………………..
Chapter 2 SATELLITE FORMATION-KEEPING USING THE
FUNDAMENTAL EQUATION IN THE PRESENCE OF ATTITUDE
REQUIREMENTS……………………………………………………..
2.1 Fundamental Equation……..………………………………………………
2.2 Formation-Keeping Equations of Motion……………………….……….
2.2.1 Formation-keeping control of the satellites in an unperturbed
circular orbit………………………………………………………...
2.2.2 Solving the problem directly in the Hill frame assuming an arbitrary
reference orbit of the leader satellite………………………………..
2.2.3 In case of incorrect initial conditions………………………………..
2.2.4 Effect of the nonuniform Earth gravity perturbation………………..
2.2.5 Formation-keeping with attitude constraints………………………..
2.3 Numerical Example……………………………………………………...
Chapter 3 SEPARATION PRINCIPLE: EXPLICIT DETERMINATION OF
THE CONTROL FORCE AND TORQUE FOR SATELLITE
FORMATION-KEEPING…………………………………………….
3.1 Determining the Explicit Form of the Control Forces and Torques…..…
3.2 Separation Principle……………………..……………………………….
3.3 Numerical Example……..………………………………………………..
Chapter 4 METHODOLOGY FOR SATELLITE FORMATION-KEEPING IN
THE PRESENCE OF SYSTEM UNCERTAINTIES……………..…
4.1 Introduction………………………………………………………………
4.2 The Description of Constrained Mechanical Systems………….………
4.3 Formation-Keeping Equations of Motion: The Controlled Nominal
System……….………………………………………………………..…
ii
vi
ix
x
1
1
4
8
8
11
16
16
19
19
26
30
33
37
41
59
59
75
81
98
98
99
104
v
4.4 Formation-Keeping Equation of Motion: The Controlled Actual System.
4.5 Generalized Sliding Surface Controller…………………………………..
4.6 Numerical Results and Simulations……………………………………...
Chapter 5 NEW SOLUTIONS TO THE INVERSE PROBLEM OF THE
CALCULUS OF VARIATIONS IN A SINGLE-DEGREE-OF-
FREEDOM SYSTEM.....................................................................…….
106
109
117
130
5.1 The Solution to the Inverse Problem………...……………………………… 130
5.2 Variational Principle for Systems with Nonpotential Forces……………. 135
Chapter 6 LAGRANGIANS FOR DAMPED LINEAR MULTI-DEGREE-OF-
FREEDOM SYSTEMS………...............................................................
6.1 Lagrangians for Damped Linear Two-Degree-of-Freedom Systems…………
6.2 Lagrangians for Special Constant Coefficient Multi-Degree-of-Freedom
Linear Systems…………………………………………………………...
Chapter 7 GENERAL CONCLUSIONS……………………….…….....................
Chapter 8 FUTURE WORK……………………………………………………….
Bibliography………………………………………………………………………..
Appendix A…………………………………………………………………………
Appendix B…………………………………………………………………………
Appendix C…………………………………………………………………………
Appendix D…………………………………………………………………………
Appendix E………………………………………………………………………….
143
143
160
172
176
178
183
186
189
192
195
vi
List of Figures
Figure 2-1: ECI frame (X-Y-Z) and Hill frame (x-y-z). The distance between the
leader satellite and the follower satellite must be maintained to be ρ in
the yz-plane of the local, rotating coordinate frame x-y-z ……………. 21
Figure 2-2: The leader satellite’s trajectory in the ECI frame…………………….. 43
Figure 2-3: Time history of the altitude and the latitude of the leader satellite…… 44
Figure 2-4: Controlled motion in the yz- and xz-planes of the first follower satellite.
The follower satellite oscillates 4 times about the circle of radius ρ
0
for
each revolution around the leader satellite (as seen in the Hill frame)... 50
Figure 2-5: Controlled motion in the yz- and xz-planes of the second follower
satellite………………………………………………………………… 50
Figure 2-6: Uncontrolled motion in the yz- and xz-planes of the first follower
satellite……………………………………………………………….... 51
Figure 2-7: Uncontrolled motion in the yz- and xz-planes of the second follower
satellite………………………………………………………………… 51
Figure 2-8: Quaternion of the first follower satellite…………………………...…. 52
Figure 2-9: Quaternion of the second follower satellite…………………………... 52
Figure 2-10: Required control forces per unit mass in the Hill frame and the force
magnitude for the first follower satellite……………………………. 55
Figure 2-11: Required control forces per unit mass in the Hill frame and the force
magnitude for the second follower satellite……………………….. 55
Figure 2-12: Required control torques per unit mass in the body frame and the torque
magnitude for the first follower satellite…………………………….. 56
Figure 2-13: Required control torques per unit mass in the body frame and the torque
magnitude for the second follower satellite………………………..... 56
Figure 2-14: Errors for the first follower satellite in the satisfaction of the six
trajectory requirements……………………………………………… 57
vii
Figure 2-15: Errors for the second follower satellite in the satisfaction of the six
trajectory requirements……………………...………………………. 57
Figure 2-16: Errors in the first follower satellite’s trajectory requirements over the
time duration 1.95 P
L
– 2 P
L
…………………..……………………. 58
Figure 2-17: Errors in the second follower satellite’s trajectory requirements over the
time duration 1.95 P
L
– 2 P
L
……………….………..……………… 58
Figure 3-1: The leader satellite’s trajectory in the ECI frame…………………….. 83
Figure 3-2: Controlled orbital motion in the yz- and xz-planes of the follower
satellite. The follower satellite stays on the circle of radius ρ around
the leader satellite (as seen in the Hill frame)………………………… 87
Figure 3-3: Uncontrolled orbital motion in the yz- and xz-planes of the follower
satellite (as seen in the Hill Frame). The initial and terminal points
indicate the start and end of the integration interval………………… 87
Figure 3-4: Quaternion of the follower satellite with control……………………... 88
Figure 3-5: Quaternion of the follower satellite without control………………….. 88
Figure 3-6: Required control forces per unit mass in the ECI frame for the follower
satellite from Eq. (3.46)……………………………………………….. 92
Figure 3-7: Required control torques per unit mass in the body frame for the follower
satellite from Eq. (3.46)……………………………………………….. 92
Figure 3-8: Errors in the satisfaction of the orbital and attitude requirements for the
follower satellite with control in the presence of sensor measurement
errors showing asymptotic convergence to the desired orbital and
attitude requirements………………………………………………….. 93
Figure 3-9: Errors in the orbital and attitude requirements when no control forces
and no control torques are applied to the follower satellite…………… 93
Figure 3-10: Each component of the P
1
A
21
T
α
O,A
to check if the separation principle
holds……………………………………………………………………. 94
Figure 3-11: Block diagram showing compensation for unmodelled dynamics using
a PD controller for robust control…………………………………… 96
Figure 3-12: Tracking errors in position showing robust control in the presence of
small unmodelled dynamics………………………………………… 96
viii
Figure 3-13: Control forces per unit mass exerted by the PD controller………….. 97
Figure 4-1: Constrained motion of the nominal system with no uncertainties
assumed……………………………………………………………….. 121
Figure 4-2: Time history of quaternions of the nominal system with no uncertainties
assumed……………………………………………………………….. 122
Figure 4-3: Required control forces to satisfy the nominal orbital constraints……. 122
Figure 4-4: Required control torques to satisfy the nominal attitude constraints…. 123
Figure 4-5: The errors in the satisfaction with the nominal constraints…………... 123
Figure 4-6: Motion of the actual system when ±10% uncertainties in the mass and
moments of inertia of the follower satellite are involved……………... 125
Figure 4-7: Quaternions of the actual system when ±10% uncertainties in the mass
and moments of inertia of the follower satellite are involved………… 125
Figure 4-8: Orbital errors between the nominal and controlled systems………….. 128
Figure 4-9: Attitude (quaternion) errors between the nominal and controlled
systems………………………………………………………………… 128
Figure 4-10: Required additional control forces to compensate for uncertainties in
the mass and moments of inertia of the follower satellite…………….. 129
Figure 4-11: Required additional control torques to compensate for uncertainties in
the mass and moments of inertia of the follower satellite…………….. 129
Figure 6-1: Linear 2-DOF mass-spring-damper system with b = 2m…………..…. 146
Figure 6-2: General 2-DOF mass-spring-damper system…………………………. 155
Figure 6-3: General 6-DOF mass-spring-damper system…………………………. 168
Figure C-1: During attitude control the angle θ remains small, which guarantess
nonzero
3
T ξ ………………………………………………………….... 191
ix
List of Tables
Table 3-1: Various conditions for the separation principle depending on the rank s
of A
11
……………………………………………………………………. 81
Table 5-1: The function , , G t x t x t
of Eq. (5.2) and the corresponding function
, , f t x r of Eq. (5.10)………………………………………………….. 140
Table 6-1: Three cases when a Lagrangian function of the system described by
Eqs. (6.8) exists and a corresponding Lagrangian for each case………... 155
x
Abstract
This dissertation handles two topics that are fundamental in understanding
aerospace and mechanical systems. In the first topic, a new controller for
maintaining a desired formation of multiple satellites is developed in the presence of
attitude tracking requirements under various perturbations and model uncertainties.
Both orbital and attitude dynamics are simultaneously considered and continuous
thrust propulsion systems are assumed. Unlike most studies, no linearizations and/or
approximations are made in dynamics or the controllers. The novel controller
developed herein provides a remarkable improvement over the current state-of-the-
art in that the control force and torque to satisfy the given orbital and/or attitude
constraints are obtained in completely closed form. This new, analytical solution can
be easily used for on-orbit, real-time control with low computational burden.
Furthermore, it is useful in estimating the magnitude of the required control inputs
and results in some interesting consequences, including the separation principle,
which describes when the orbital control can be completely decoupled from the
attitudinal control required to be applied to the follower satellite. Also, an additional
additive controller that compensates for model uncertainties is developed. This is
done by using the desired trajectory of the nominal system as the tracking signal,
xi
and is based on a generalization of the concept of sliding surfaces. The resulting
closed-form control causes the desired attitude and orbital requirements of the
nominal system to be met in the presence of unknown, but bounded, model
uncertainties.
The second topic deals with the inverse problem of the calculus variations. To
develop the equations of motion for a mechanical system, it is usual to define a
Lagrangian function first and then obtain the equations of motion by applying the
Lagrange equation when the forces are derivable from a potential function.
However, it is interesting and useful to investigate the types of forces that can
engendered when substituting the forces into the Lagrange equation guarantees the
proper equations of motion even when the forces do not arise from a potential, that
is, when they are nonconservative. By extending the existing Bolza’s approach, a
more general and systematic way of handling the problem is found for single-
degree-of-freedom systems. Then, new results for linearly damped multi-degree-of-
freedom systems are obtained using extensions of results for single-degree-of-
freedom systems. Conservation laws for these damped multi-degree-of-freedom
systems are also found using the Lagrangians obtained. These new findings can be
used for development of approximate solutions of nonlinear differential equations or
various numerical techniques.
1
Chapter 1
GENERAL INTRODUCTION
1.1 Motivation
This dissertation deals with two topics that have attracted significant interest in recent
years. The first topic deals with the problem of precision formation flying of satellites
including precision attitude control--the so-called full body problem in dynamics--, and
the second topic deals with the inverse problem of the calculus of variations.
A. Over the past few years there has been a growing interest in satellite formation
flying. Formation flying has been extensively investigated as a means to expand
the capabilities of existing space missions. Technologies, like interferometry,
associated with satellites flying in a formation demand precise formation-keeping
and hence this area of research has grown enormously in the last decade or so. To
successfully achieve a mission goal it is often essential to accurately maintain a
formation geometry as desired even under a variety of perturbations. Formation
flying is usually described in terms of a real or fictitious “leader” satellite that
serves as a reference to the others in the formation and a set of “follower”
satellites that maintain a prescribed configuration around the leader satellite while
moving in some fixed manner relative to it. Previous work mainly focuses on
2
linearized approaches to the formation-keeping problem in which the trajectory of
the leader satellite is taken to be a simple circular or elliptical Keplerian orbit.
When nonlinear dynamics is considered, the controllers for formation-keeping
have been obtained only via the Lyapunov method or by numerical approaches.
Also, only a few researchers incorporated attitude dynamics in the formation-
keeping problem. In a real formation flying mission, it is important for a satellite
to keep a specific orientation for the sake of stable power supply, successful
mission goal achievement, and so on. Furthermore, the designed controller
assumes the exact measurement of the states in real time, but it is very hard or
impossible to exactly estimate them due to unmodeled dynamics, sensor noise,
and so on. In this dissertation, hence, a simple analytical approach for formation-
keeping in the presence of attitude requirements under model uncertainties is
presented. More specifically, assuming a variety of perturbations in space, a
nonlinear exact controller is developed in a closed form even when the leader
satellite follows an arbitrary reference orbit revolving around the nonspherical
Earth, and the follower satellites are to point towards a specific spot in space
during the maneuver. Also, initial insertion errors are considered and it is shown
that these errors are made to asymptotically converge to zero as time progresses
when the newly obtained controller developed in this dissertation is implemented.
3
B. Regarding the second topic, to develop the equations of motion for a mechanical
system from Hamilton’s principle, it is usual to define a Lagrangian function from
the kinetic and potential energies first and then obtain the equations of motion by
applying the Lagrange equation when the constraints on the system are holonomic
and the forces are derivable from a potential function. However, it is interesting to
investigate the types of forces that can be engendered when substituting the forces
into the Lagrange equation guarantees the proper equations of motion even when
the forces acting on a system do not arise from a potential, that is, when they are
nonconservative. This is often called the inverse problem of the calculus of
variations or the inverse problem for Lagrangian dynamics. Many researchers
have been attracted to this problem that has served as a basis for development of
approximate methods of analysis and various numerical techniques, but previous
work has mainly focused on a single degree-of-freedom (SDOF) system and
necessitates experience, intuition, or both. In the current dissertation, some more
general nonpotential systems are included, and a more systematic way of handling
the inverse problem of the calculus of variations is provided. Also, new results for
linear multi-degree-of-freedom (MDOF) systems are provided using extensions of
results for SDOF systems.
4
1.2 Background
For formation-keeping problem of satellites, traditionally linear control theories have
been used, based on linearized equations of relative motion such as a Hill-Clohessy-
Wiltshire equation (Hill [1], Clohessy and Wiltshire [2]) for a circular reference orbit or a
Tschauner-Hempel equation (Tschauner and Hempel [3]) for an elliptical reference orbit,
or the J
2
-perturbed linearized equation (Schweighart and Sedwick [4]). Yan et al. [5]
designed a pulse-based controller for the satellites’ periodic motion using the linear
quadratic regulator (LQR) technique. Sparks [6] studied the same problem and
considered nonspherical gravity perturbations. Vaddi and Vadali [7] presented Lyapunov,
LQR, and period-matching control schemes that generate a projected circular orbit
(PCO), which is also chosen as the desired formation geometry in this dissertation. Since
the solutions based on the linearized equations about Keplerian orbits are suitable only
for very small-sized formations due to the linearizations used, Milam et al. [8] used a
numerical solver, called the Nonlinear Trajectory Generation software package to find an
optimal controller for formation-keeping. Based on nonlinear dynamics describing the
relative motion around a circular reference orbit, De Queiroz et al. [9] developed an
adaptive controller for tracking the relative position.
A few researchers incorporated attitude dynamics in formation flying. Kristiansen et
al. [10] derived a nonlinear model describing relative translational and rotational motion.
5
Long and Hall [11] designed a controller that uses Lyapunov control theory for target
tracking while countering the effect of gravity gradient torques. Ahn and Kim [12]
assumed a formation of satellites as a virtual rigid-body structure to develop an algorithm
for pointing to a specified target. Assuming that the satellites are in unperturbed
Keplerian orbits, Segal and Gurfil [13] derived more general description of the relative
translational motion between two arbitrary points which may not coincide with the
satellites’ centers of mass.
Previous work is thus seen to mainly focus on linearized approaches in which the
trajectory of the leader satellite is taken to be a simpler circular or elliptical Keplerian
orbit. When nonlinear dynamics is considered, the controllers for formation-keeping have
been obtained via the Lyapunov method or by numerical approaches. Recently, a
nonlinear controller that preserves the nonlinearities of the governing dynamics was
proposed by Udwadia and Kalaba (Udwadia and Kalaba [14-23]). Due to its simplicity
and generality, this method has been applied to highly constrained problems such as
control of tethered satellites (Schutte and Dooley [24]) or orbit transfers under
perturbations (Alizadeh and Villac [25]). For the formation-keeping problem, Udwadia et
al. [26] solved the precision tumbling and precision tracking problem for a nonlinear,
nonautonomous multi-body system orbiting a central body with a non-uniform gravity
field. They included attitude dynamics and used a closed-form control methodology to
6
achieve exact control. The inertial frame was used and the vibrating nonlinear multi-body
system was required to be in a precise circular orbit around the central body, while it
underwent arbitrarily prescribed on-orbit precision tumbling. But this paper (Udwadia et
al. [26]) used a simple, circular reference orbit and no initial insertion errors were
considered. A few researchers considered parameter uncertainties or assumed a lack of
exact measurement information. [27-29] However, most of them did not consider any
attitude dynamics and the resulting steady-state errors are not still small enough.
For the inverse problem of the calculus of variations, Bolza [30] gave a general
procedure for finding a Lagrangian for a SDOF system. This was followed by Douglas
[31,32] who obtained the necessary and sufficient conditions for the existence of a
Lagrangian for MDOF systems. However, it is difficult to obtain the Lagrangians for
given, specific forces. Leitmann [33] provided some examples of nonpotential forces and
the corresponding Lagrangians for which a variational principle exists. A SDOF system
was considered and two examples were provided. Recently, the so-called semi-inverse
method (He [34]) has attracted much attention due to its simplicity and applicability to
certain cases. However, this method is not general in that it assumes a specific form for
the Lagrangian function which needs to be obtained by chance or from intuition.
However, [30], [33], and [34] consider only SDOF systems, and the analysis of this case
is relatively easy because Darboux [35] proved in the nineteenth century that a
7
Lagrangian can be always found for the inverse problem for such SDOF systems. For
MDOF systems the configuration variables are coupled with one another and this makes
it difficult to solve the inverse problem. The general conditions for the existence of
Lagrangians were apparently first obtained by Helmholtz (Helmholtz [36], Santilli [37])
and are usually referred to as the Helmholtz conditions. Later, Douglas [31,32] analyzed
in great detail the case of two degrees of freedom (2-DOF) and obtained the necessary
and sufficient conditions for their existence without utilizing these conditions. Using
Douglas’ results Hojman and Ramos [38] proposed a simpler method to determine the
existence of a Lagrangian for two-dimensional problems in which the potential function
does not explicitly contain the generalized velocities. Mestdag et al. [39] derived the
conditions under which there exists a Lagrangian and a dissipation function on the right
hand side of the more general form of the Euler-Lagrange equation. They also provided
some nonconservative systems to which their approach can be applied.
Although they are explicit and provide the exact conditions for the existence of a
Lagrangian for a MDOF system, the Helmholtz conditions [36,37] are near-impossible to
solve for general n-degree-of-freedom (n-DOF) systems (and most often even for 2-
degree-of-freedom (2-DOF) systems), and from a practical standpoint they provide little
assistance in solving the inverse problem for such systems. Instead, it is sometimes
desirable to develop a simpler method that circumvents the use of the Helmholtz
8
conditions and that is easily applicable to MDOF systems commonly arising in the theory
of vibrations.
1.3 Research Objectives
The two major challenges for this dissertation are:
(1) To explicitly determine the control forces and torques for satellite formation-
keeping with attitude requirements under various perturbations in space and to
develop an additional additive controller that compensates for model
uncertainties.
(2) To provide a systematic way to solve the inverse problem of the calculus of
variations and to obtain the Lagrangians for MDOF systems in a simple manner,
using insights obtained from the understanding of the inverse problem for SDOF
systems.
1.4 The Study’s Contributions
This dissertaion is divided into two main sections. The purpose of the first section is to
present exact formation-keeping guidance schemes that determine control inputs and state
trajectories. The solution presented in this dissertation includes attitude constraints and is
new in that 1) we assume that the leader satellite’s reference orbit can be any arbitrary
one, 2) a method to correct initial-condition misalignment problems is provided, 3) the
9
new solution leads to some interesting consequences, including the separation principle
that points out the circumstances under which the orbital control is completely unaffected
by the attitude control, and 4) the effect of model uncertainties is also considered and it is
shown by examples that the nonlinear controller proposed herein when augmented by an
additional additive controller based on a generalization of the notion of sliding surface
control can satisfy the control requirements with considerable accuracy in the presence of
such uncertainties. New solutions are obtained that consider all of these complex
situations, and it seems that no previous research has accomplished this.
First, an arbitrary leader satellite’s orbit is assumed and no restrictions on the orbit are
placed. That is, the reference orbit can be an arbitrary function of time to include the
effects of various disturbances and/or thrust. Next, to extend the practical utility of the
exact, explicit closed-form solutions obtained in this dissertation, a method to correct
initial-condition misalignment problems is included, since in practice there is no
guarantee that the follower satellites are initially inserted exactly into the desired orbits.
Also, it is investigated under what conditions the attitude dynamics can be ignored, or
separated away, when designing a controller for exact orbital control, leading to the
statement of the Separation Principle. As a special case, if there are three independent
orbital constraints given, the separation principle is always satisfied throughout the
maneuver. Lastly, an additional additive controller that compensates for model
10
uncertainties is developed. This is done by using the desired trajectory of the nominal
system as the tracking signal, and is based on a generalization of the concept of sliding
surface control. This methodology is verified by an example in which a follower satellite
has uncertain mass and moments of inertia, and numerical results show high accuracy of
the control methodology in maintaining the desired formation requirements even when
model uncertainties are present in the system.
In this dissertation, both orbital and attitude dynamics are compactly handled in a
unified way, and the approach shows the possibility of precisely controlling highly
nonlinear systems. A numerical example shows that the controller designed in the current
dissertation guarantees successful tracking of the given constraint trajectories with high
precision under a variety of pertuations and model uncertainties.
In the second section, the results of the previous research (Leitmann [33]) are
extended to some more general nonpotential systems, and a more systematic way of
handling the inverse problem of the calculus of variations in SDOF systems is provided.
The examples given in Leitmann [33] arise as special cases of the results provided herein.
Then, the general results are applied to some specific systems to indicate the nature of the
nonpotential forces which the results encompass. Next, new results for linearly damped
MDOF systems are obtained using extensions of results for SDOF systems. The solution
to the inverse problem for an n-DOF linear gyroscopic system is obtained as a special
11
case. MDOF systems that commonly arise in linear vibration theory with symmetric
mass, damping, and stiffness matrices are handled similarly in a simple manner.
Conservation laws for these damped MDOF systems are found using the Lagrangians
obtained, and several examples are provided. These new findings can be used for
development of approximate solutions of nonlinear differential equations or various
numerical techniques.
1.5 Organization
This study is organized into the following eight chapters:
Chapter 2: Satellite Formation-Keeping Using the Fundamental Equation in the
Presence of Attitude Requirements
In this chapter, a simple analytical method for the formation-keeping problem is
presented. At first, the simplest case where an unperturbed, circular reference orbit of the
leader satellite is considered. Also, the problem is solved in the inertial frame first, then
transformed into the Hill frame. This step is a stepping stone for further, more general,
developments. Next, it is shown that the problem can be solved alternatively in the Hill
frame considering an arbitrary reference orbit and insertion errors. Finally, a nonuniform
gravity perturbation of the nonspherical Earth is considered, resulting in a highly
nonlinear, nonautonomous control problem. This precision control problem is handled in
12
a closed form, a real-time framework, so that autonomous control of such formation
flying can be brought within the realm of practicality. Also, one interesting fact about the
number of constraint equations is provided. For the attitude target-pointing constraint, we
should use two equations. It seems like that we could have used only one constraint
equation (instead of two) by combining the existing two constraint equations, but it is
shown that this single constraint equation does not result in correct answers and why the
single constraint equation fails is explained as well. Numerical simulations demonstrating
the efficacy and the accuracy of the analytic solutions developed in this dissertation are
also provided.
Chapter 3: Separation Principle: Explicit Determination of the Control Force and
Torque for Satellite Formation-Keeping
In this chapter, it is shown that the closed-form solutions obtained in Chapter 2 permit
some important physical insights into the nature of the control, something that would be
difficult to obtain from purely computational approaches. For example, the conditions
under which orbital control is independent of the attitude of the satellite, so that one need
not consider any attitude in designing the orbital controller, are investigated. As a special
case, if there are three independent orbital constraints given, the separation principle is
always satisfied throughout the maneuver. Finally, the effect of unmodelled dynamics is
considered and it is shown that by augmenting the closed-form nonlinear controller
13
obtained herein with a simple PD controller the control requirements can be met with
high accuracy, thereby pointing out the robustness with which unmodelled dynamics can
be easily handled.
Chapter 4: Methodology for Satellite Formation-Keeping in the Presence of System
Uncertainties
This chapter proposes a formation-keeping control methodology that includes both
attitude and orbital control requirements in the presence of model uncertainties. First, a
nominal system model that provides our best assessment of the real-life uncertain system
is defined, and a nonlinear controller is developed following the procedures provided in
Chapters 2 and 3. Next, since the nominal system assumes no uncertainties, an additional
additive controller that compensates for model uncertainties is developed. This is done by
using the desired trajectory of the nominal system as the tracking signal, and is based on
a generalization of the concept of sliding surface control. The resulting control input
causes the desired attitude and orbital requirements of the nominal system to be met in
the presence of unknown, but bounded, model uncertainties.
14
Chapter 5: New Solutions to the Inverse Problem of the Calculus of Variations in a
Single Degree-of-Freedom System
This chapter deals with the inverse problem of finding a suitable integrand so that upon
the use of the calculus of variations, one obtains the equations of motion for SDOF
systems in which the forces are nonpotential. First, a simple and systematic way to get a
nonpotential force and its corresponding Lagrangian for a SDOF system is found. The
main difficulties lie in performing the necessary integrations explicitly. Next, new
Lagrangians for some specific systems are proposed for a SDOF system. Finally, we
apply the obtained results to some specific SDOF systems to indicate the nature of the
nonpotential forces, which the results encompass.
Chapter 6: Lagrangians for Damped Linear Multi-Degree-of-Freedom Systems
In this chapter the problem of finding Lagrangians for damped, linear MDOF systems is
dealt with. New results for such MDOF systems are obtained using extensions of results
for SDOF systems which are given in Chapter 5. The solution to the inverse problem for
an n-DOF linear gyroscopic system is obtained as a special case. MDOF systems that
commonly arise in linear vibration theory with symmetric mass, damping, and stiffness
matrices are handled similarly in a simple manner. Conservation laws for these damped
15
MDOF systems are also found using the Lagrangians obtained, and several examples are
provided.
Chapter 7: General Conclusions
Chapter 7 outlines the analytic control scheme for precise formation-keeping from
Chapter 2, the separation principle from Chapter 3, and the new control methodology to
cope with model uncertainties from Chapter 4. Also, it summarizes the systematic
procedure in getting new nonpotential forces of SDOF system which are integrand
functions in the variational principle from Chapter 5 and the straightforward extension to
MDOF systems as well as the usefulness of the new method from Chapter 6.
Chapter 8: Future Work
Lastly, Chapter 8 proposes plans for future work and suggests follow-up studies that
could be based on this dissertation.
16
Chapter 2
SATELLITE FORMATION-KEEPING USING THE
FUNDAMENTAL EQUATION IN THE PRESENCE OF
ATTITUDE REQUIREMENTS
2.1 Fundamental Equation
This section deals with the explicit equation of motion for a constrained mechanical
system. First, we consider an unconstrained, discrete dynamical system of N rigid bodies
(or satellites in this dissertation) in generalized coordinates. The generalized
displacement vector and velocity vector are denoted by t q and t q respectively,
where
1 2
,
n
t q t q t q t
T
q
1 2
,
n
t q t q t q t
T
q
(2.1)
t represents time, the superscript “ T ” denotes the transpose of a vector or a matrix, and
n is the number of the generalized coordinates. If the given forces impressed on the
bodies are denoted by
1 2
,
n
t Q t Q t Q t
T
Q
(2.2)
then the Lagrange equation for the unconstrained motion of the system at a certain instant
of time t can now be expressed as an n by 1 matrix equation
, , t t t t Mq Q q q
(2.3)
17
or
, , : , t t t t t
-1
q M Q q q a
(2.4)
where M (which is positive definite in this dissertation) is an n by n symmetric mass
matrix and t a is the n by 1 unconstrained acceleration (vector) of the system at the
time t .
It is now considered that the unconstrained system described by Eqs. (2.3) and (2.4) is
subjected to p constraints which are of the form
, , 0, 1,2, , .
j
t j p q q (2.5)
These p constraints include both orbital and attitude requirements for precise formation-
keeping in the current dissertation. Differentiating Eq. (2.5) with respect to time once (for
nonholonomic constraints) or twice (for holonomic constraints), we get the following
constraint equation
, , , , , t t t t t t A q q q b q q
(2.6)
where the matrix A is a p by n matrix and b is a p by 1 vector.
The presence of the constraints, Eq. (2.6), causes additional forces t
C
Q to be
applied to the bodies and the resulting equation of motion becomes
, , , t t t t t
C
Mq Q q q Q
(2.7)
where the (control) force t
C
Q makes the system exactly satisfy the given constraints of
Eq. (2.6) at each instant of time. Historically there have been a number of attempts to
18
determine t
C
Q , and recently a simple and insightful approach was proposed by
Udwadia and Kalaba (Udwadia and Kalaba [14]), which gives the constraint force t
C
Q
explicitly as
, t
+ +
C 1/2 -1/2 T -1 T
Q M AM b Aa A AM A b Aa (2.8)
where the superscript “+” represents the Moore-Penrose generalized inverse.
Thus, the explicit equation of motion of the constrained system is the following
+ +
1/2 -1/2 T -1 T
Mq Ma M AM b Aa Ma A AM A b Aa (2.9)
or
.
+ +
-1/2 -1/2 -1 T -1 T
q a M AM b Aa a M A AM A b Aa (2.10)
In this dissertation, Eq. (2.10) is referred to as the fundamental equation, and the
constraint force t
C
Q in Eq. (2.8) is used for controlling the follower satellites so that
each of them satisfies its formation-keeping requirements. Furthermore, it is shown in
Udwadia [20] that this constraint force t
C
Q minimizes at each instant of time t the
quantity , t t t
T
C -1 C
Q M q Q
.
In what follows we shall use the term “unconstrained motion” to mean “uncontrolled
motion”; the term “constrained motion” to mean “controlled motion”; and the term
“constraint force” to mean “control force.”
19
2.2 Formation-Keeping Equations of Motion
It is assumed in this dissertation that there are N follower satellites in a formation where
either a real or a fictitious leader satellite is located at its center. The i
th
follower satellite
has a mass
( ) i
m and has a diagonal inertia matrix,
( ) i
J , in which the moments of inertia
along its body-fixed principal axes of inertia are placed.
First, let us start with the simplest case.
2.2.1 Formation-keeping control of the satellites in an unperturbed circular orbit
We assume that the leader satellite’s orbit is circular. Also, there is no perturbation
applied on the satellites, and the given constraints are satisfied from the initial time. The
inertial orbital motion of the follower orbiting the spherical Earth is governed by the
relation (Vallado [40])
3/2
2 2 2
,
ECI
X X
GM
Y Y
X Y Z
Z Z
a
(2.11)
where X Y Z
T
is the position vector of the center of mass of the follower satellite in
the ECI frame, G is the universal gravitational constant, and M
is the mass of the
Earth. The subscript “ECI” denotes that Eq. (2.11) is described in the ECI frame. This
inertial coordinate frame originates at the center of the Earth, the X axis points towards
20
the vernal equinox, the Y axis is 90
to the east in the equatorial plane, and the Z axis
extends through the North Pole. (Fig. 2-1)
Next, we specify two constraints on the follower satellite’s trajectory. As a first step,
it is assumed that there is only one follower satellite in a formation. Our goal is to
maintain the distance between the leader satellite and follower satellite to be constant, say
, when projected on the local horizontal plane. Also, to prevent the x coordinate from
diverging, another constraint between x and z coordinates is applied. To describe these
constraints, it is convenient to define a new coordinate frame, called the Hill frame or the
LVLH (Local-Vertical and Local-Horizontal) frame. The origin of this frame is located at
the leader, the x -axis lies in the radial direction, the y -axis is in the along-track
direction, and the z -axis lies along the orbital angular momentum vector, thus
completing a right-handed coordinate system. (See Fig. 2-1)
21
Fig. 2-1 ECI frame (X-Y-Z) and Hill frame (x-y-z). The distance between the leader satellite and the
follower satellite must be maintained to be ρ in the yz-plane of the local, rotating coordinate frame x-
y-z.
In this Hill frame, let us choose the following two constraints:
2 2 2
, y z (2.12a)
2 0, x z (2.12b)
where is the constant distance between the leader satellite and the follower satellite.
Eq. (2.12b) is added to make the relative motion to be bounded in every axis. Since Eqs.
(2.12) are holonomic constraints, we differentiate Eqs. (2.12) twice with respect to time
and the equations become:
22
2 2
, yy zz y z (2.13a)
2 0, x z (2.13b)
where x , y , and z are described in the Hill frame. It is noted that Eqs. (2.13) are of the
form Eq. (2.6) where
2 2
0
, , .
2 0 1 0
x
y z y z
y
z
A q b
(2.14)
Since the Hill frame is a rotating or noninertial frame, we transform Eq. (2.13) into the
one described in the ECI (inertial) frame. But it will be shortly shown that the problem
can be directly solved in the Hill frame. From the definition of the Hill frame, it is easy to
show that the transformation equation which maps the ECI frame ( , , ) X Y Z to the Hill
frame ( , , ) x y z is given by:
.
L
x r X
y Y
z Z
= R
(2.15)
Here,
L
r is the radius of the leader satellite’s circular orbit and the transformation matrix
R is the following:
cos( )cos( ) sin( ) cos( )sin( ) sin( )cos( ) cos( ) cos( )sin( ) sin( )sin( )
cos( )sin( ) sin( )cos( )cos( ) cos( ) cos( ) cos( ) sin( )sin( ) sin( ) cos( ) ,
sin( )sin( ) cos( )sin( ) cos( )
i i i
i i i
i i i
R
(2.16)
where , i , and are the right ascension of the ascending node, the inclination, and the
argument of perigee of the leader satellite, respectively. The leader satellite’s orbit is
circular, i.e. the perigee is not defined, so it follows that is used as the angle between
23
the ascending node and the leader satellite’s position vector. Since the leader satellite is
assumed to be in an unperturbed circular orbit, and i are constant with time, and is
a linear function of time. After differentiating Eq. (2.15) twice with respect to time and
substituting into Eqs. (2.13) (the constraint equations in the Hill frame), the following
constraint equation in the ECI frame is obtained:
11 12 13 1
21 22 23 2
.
X
A A A b
Y
A A A b
Z
(2.17)
Then, each component of the matrices A and b can be written as follows:
2
2 2
11
2
cos( )sin( ) sin( )cos( )cos( ) sin ( )sin ( )
cos( )sin( ) sin( )cos( )cos( ) sin( )sin( ) cos( )cos( )cos( ) sin( )cos( )sin ( )
sin( )sin( )cos( ) sin( )cos( ) cos( )sin( ) sin( )co
A X i i
Y i i i
Z i i i
s( )cos( ) , i
(2.18a)
2
12
2
2 2
cos( )sin( ) sin( )cos( )cos( ) sin( )sin( ) cos( )cos( )cos( ) sin( )cos( )sin ( )
cos( )cos( )cos( ) sin( )sin( ) cos ( )sin ( )
sin( )cos( ) cos( )cos( )cos( ) sin( )sin( ) cos( )si
A X i i i
Y i i
Z i i
n( )cos( ) , i i
(2.18b)
13
2 2 2
sin( )sin( )cos( ) sin( )cos( ) cos( )sin( ) sin( )cos( )cos( )
sin( )cos( ) cos( )cos( )cos( ) sin( )sin( ) cos( )sin( )cos( )
sin ( )cos ( ) cos ( ) ,
A X i i i i
Y i i i i
Z i i
(2.18c)
and
21
2cos( )cos( ) 2sin( )cos( )sin( ) sin( )sin( ), A i i (2.19a)
22
2sin( )cos( ) 2cos( )cos( )sin( ) cos( )sin( ), A i i (2.19b)
23
2sin( )sin( ) cos( ). A i i (2.19c)
24
In addition,
1
- - 2
- - 2 - ,
T T
T T T
X X
b X Y Z Y X Y Z Y
Z Z
X X X
X Y Z Y X Y Z Y X Y Z Y
Z Z Z
R R R R
R R R R R R
(2.20a)
and
2
2 1 4 2 ,
X X
b Y Y
Z Z
R R
(2.20b)
where R is a 2 by 3 matrix that is composed with the second and third rows of the
matrix R (see Eq. (2.16)) and R
is a 2 by 3 matrix that is composed with the first and
third rows of the matrix R. The generalized Moore-Penrose inverse of the matrix A
can be obtained in the following closed form:
( ) ,
+ + + +
1
A A c c (2.21)
where
11
12 2 2 2
11 12 13
13
1
,
A
A
A A A
A
+
1
A =
(2.22)
11 21 12 22 13 23
2 2 2
11 12 13
,
A A A A A A
A A A
(2.23)
25
21 11
22 12 2 2 2 2 2 2 2
21 22 23 11 12 13
23 13
1
.
( )
A A
A A
A A A A A A
A A
+
c (2.24)
Finally, substituting all these terms into the fundamental equation (Eq. (2.10)), we
explicitly obtain the equation of motion for the follower:
1
3/2
2 2 2
2
11 12 13
3/2
2 2 2
21 22 23
( )
( ) .
X
b
GM
Y
b
X Y Z
Z
X
A A A
GM
Y
A A A
X Y Z
Z
+
-1/2 -1/2
+ + +
1
+ + +
1
x = a M AM b Aa
= A c c
A c c
(2.25)
The control force that is needed to satisfy the constraints can be explicitly obtained:
11 12 13 1
3/2
2 2 2
21 22 23 2
( ) ( ) ,
X
A A A b GMm
m Y
A A A b
X Y Z
Z
+
c 1/2 -1/2
+ + + + + +
1 1
Q = M AM b Aa
= A c c A c c
(2.26)
where m is the mass of the follower satellite. In Eqs. (2.25) and (2.26),
+
1
A , , and
+
c
are explicitly given in Eqs. (2.22)-(2.24) in closed form. Also, the acceleration in Eq.
(2.25) and the control force in Eq. (2.26) are described in the ECI frame.
26
2.2.2 Solving the problem directly in the Hill frame assuming an arbitrary
reference orbit of the leader satellite
In the previous subsection, we solved the formation keeping problem in the inertial (ECI)
frame. However, when it comes to formation flying missions, it is sometimes preferred to
solve the problem directly in the Hill frame because the calculation becomes much more
concise and we can directly use the measured data. In this subsection, our aim is to
handle the problem wholly in the Hill frame. Let us start with the transformation
relationship between the ECI frame and the Hill frame described by Eq. (2.15):
,
L
x r X
y Y
z Z
= R
(2.15)
or
,
L L L
X x r x r x r
Y y y y
Z z z z
-1 T
R R S
(2.27)
where x y z
T
denotes the coordinates of an arbitrary point in the Hill coordinate
frame, X Y Z
T
denotes its coordinates in the ECI coordinate frame, R is an
orthogonal rotation matrix, S is the transpose of R , and
L
r (the subscript “L” denotes
quantities pertinent to the leader satellite) is the distance between the origin of the ECI
frame and the origin of the Hill frame (leader satellite). If the leader satellite’s orbit is an
unperturbed Keplerian orbit as in the previous subsection, the rotation matrix R is
27
obtained in an easy manner. If the orbit is an arbitrary function of time, however, things
are more complicated. Let us assume that the position and velocity vector of the leader
satellite are given by:
( ), ( ), ( ) , ( ), ( ), ( ) .
L L L L L L L L
X t Y t Z t X t Y t Z t
T
T
r v
(2.28)
It should be noted that these vectors are represented in the ECI frame. Then, the
instantaneous orbital angular momentum vector
L
h of the leader is given by
.
L L L
h r v (2.29)
Now, we find the rotation matrix R that maps the ECI frame to the Hill frame, that
is,
.
L
x r X
y Y
z Z
R
(2.30)
To get the matrix R , let us recall that we have three vector equations
1 0 0
ˆ ˆ
ˆ ˆ 0 , 1 , 0 ,
0 0 1
L L L L
Rr R h r Rh
(2.31)
where ˆ
L
L
L
r
r
r
and
ˆ
L
L
L
h
h
h
(unit vectors). Then, we have
1 0 0
ˆ ˆ
ˆ ˆ 0 1 0 .
0 0 1
L L L L
R r h r h
(2.32)
Since R is orthogonal, from Eq. (2.32) we obtain
28
ˆ
ˆ ˆ ˆ
ˆ ˆ ˆ .
ˆ
L
L L L L L L
L
T
T T
T
r
R r h r h h r
h
(2.33)
Each element of the matrices R , R
, and R
is given in the Appendix A, which will be
used for computing the unconstrained (uncontrolled) acceleration a in the Hill frame.
Next, let us obtain the uncontrolled (unconstrained) acceleration a (Eq. (2.11)) in the
Hill frame. Eq. (2.27) can be written more succinctly as follows:
, X Sx Sγ Sx w (2.34)
where x y z
T
x , X Y Z
T
X , 0 0
L
r
T
γ , and w Sγ . To get the
uncontrolled (unconstrained) acceleration
Hill
a in the Hill frame, we start from the
Lagrange equation
, , ,
d T T
t
dt
Q x x
x x
(2.35)
where T is the kinetic energy of the system and , ,t Q x x denotes the generalized force
of the system which includes the potential energy. Since the kinetic energy T must be
represented in the inertial frame, it follows that
1
2
1
2
1
2
1
,
2
T
T
T
T T T
T T T T T
X MX
wMSx wMSx w Mw
x SMSx x SMSx x SMw
xS MSx x S MSx x S Mw
(2.36)
29
where Eq. (2.34) has been used. Inserting Eq. (2.36) into Eq. (2.35), we get the following
equation of motion in the Hill frame
2 , , , t
T T T T T
S MSx S MSx S MSx S Mw S F X X
(2.37)
where
, ,t F X X
is given by
3/2
, , .
GM
t
T
F X X F X MX
X X
(2.38)
From Eq. (2.37), the matrix
T
S MS serves as the mass matrix in the Hill frame. However,
the matrix M is m
3 3
M I
for each particle (or satellite) where m is the mass of the
particle (or satellite), so
T
S MS is just equal to M owing to the orthogonality of the
matrix S . Finally, the uncontrolled acceleration vector a in the Hill frame can be
described as
3/2
2 ,
Hill
GM
T T T T
T
a S Sx S Sx S w x S w
Sx w Sx w
(2.39)
or more explicitly as
3/2
2
2 2
0 2 .
0
L L L L
Hill
L
x r x r x r x r
GM
y y y y
x r y z
z z z z
a RS RS
(2.40)
Having explicitly obtained the rotation matrix R given by Eq. (2.33), the uncontrolled
acceleration
Hill
a
in the Hill frame given by Eq. (2.40) can be found, and we can directly
use the explicit control force for precise formation-keeping given by using Eq. (2.8), and
the equation of motion of the controlled system by using Eq. (2.10). It must be stressed
30
here that now we can directly use the constraint equation described in the Hill frame (for
example, Eq. (2.14)), that is, we do not have to convert it to the one in the ECI frame
such as Eqs. (2.17)-(2.20), which is usually a time-consuming and troublesome task.
For instance, when the constraints Eqs. (2.13) are given, the explicit equation of
motion is obtained by
2 2
2 2
2 2
2
2 2 2
0
1
5 ,
2 0 1 5 4 0
4
z y z
y z y z
y yz
y z
z y
+
-1/2 -1/2
+
x = a M AM b Aa
= a + A b Aa
= a + a
(2.41)
where a is computed by Eq. (2.40) and
3 3
m
M I is used. Also, the control force is
explicitly given by
2 2
2 2
2 2
2
2 2 2
0
5 ,
2 0 1 5 4 0
4
m
z y z
y z y z m
y yz
y z
z y
+
c 1/2 -1/2
+
Q = M AM b Aa
= A b Aa
= a
(2.42)
which is described in the Hill frame. See the simplicity of Eqs. (2.41) and (2.42) when
compared with Eqs. (2.25) and (2.26) which has to employ the ECI frame to proceed.
2.2.3 In case of incorrect initial conditions
The constraints must be satisfied at each instant of time during the maneuver including
the initial time ( 0 t ). In practice, however, it is usually quite difficult to meet these
31
constraints at the initial time because this requires inserting the follower satellites into
orbit with the exact, required initial conditions. Hence, we need to modify our
formulation; we do this by using Baumgarte’s constraint stabilization method (Baumgarte
[41]).
In general the constraint equations may be thought of as a set of p constraints of the
form
, , 0, 1,2, , ,
i
t i p q q (2.43)
where q is an n by 1 generalized coordinates vector, n is the number of independent
variables, and t represents time. For example, in the previous subsections we have used
two (that is, 2 p ) constraints given by Eqs. (2.12):
2 2 2
1
0, y z (2.44)
2
2 0. x z (2.45)
More compactly, we can write Eq. (2.43) as follows
1
2
.
p
Φ 0
(2.46)
By differentiating the above equation with respect to time appropriately, we get the
general constraint equations mentioned earlier
, , , , . t t t A q q q b q q
(2.47)
32
Now, we consider general initial conditions that do not satisfy the given constraints,
which means (0) Φ 0 at the initial time. We modify the constraint Eq. (2.46) to
(Baumgarte [41])
, Φ αΦ βΦ 0
(2.48)
where
1 2 1 2
, , , , , , , .
p p
diag diag α β (2.49)
It is well-known that if each
i
,
i
> 0, 1,2, , i p , Φ approaches to zero
asymptotically. Thus, from Eq. (2.48), we get a more general constraint equation, which
nonetheless retains the form of Eq. (2.6).
In the previous sections, for example, we have, 2 p , so both α and β are 2 2
diagonal matrices and Eq. (2.48) yields the following constraint equation
2 2 2 2 2 2
2
1 1
0
, 2
2 0 1
2 2
x
y z yy zz y z y z
y
x z x z
z
(2.50)
which is of the form of Aq b
or Eq. (2.6). If the initial conditions do not meet the
constraints, we must utilize Eq. (2.48) instead of Eq. (2.46) and use these new quantities
A and b in Eqs. (2.8) and (2.10) respectively to get the explicit closed-form expressions
for the required control force and the explicit equation of motion of the controlled
system.
33
2.2.4 Effect of the nonuniform Earth gravity perturbation
Up to now we have assumed that there is no perturbation applied on both the leader
satellite and the follower satellites. However, this is not the case in real space missions
because satellites are exposed to various kinds of perturbation. Among these, the most
significant one applied on a satellite in low Earth orbit is the nonuniform Earth gravity
perturbation resulted from the nonspherical shape of the Earth. This perturbation makes
the problem nonautonomous and more complex, and it should be considered for more
advanced and precise control. In this subsection this perturbation will be included for
more practicality.
The inertial orbital motion of the follower satellite orbiting the nonspherical Earth is
governed by the relation (Vallado [40])
,
, 3/2
2 2 2
,
,
X ECI
ECI Y ECI
Z ECI
X X d
GM
Y Y d
X Y Z
Z Z d
a
(2.51)
where X Y Z
T
is the position vector of the center of mass of the follower in the ECI
frame, G is the universal gravitational constant, M
is the mass of the Earth, and
, X ECI
d ,
, Y ECI
d , and
, Z ECI
d are the nonuniform gravity accelerations along each of the three
(inertial) directions in the ECI frame. The subscript “ECI” is used to stress that Eq. (2.51)
is described in the ECI frame.
34
The unconstrained acceleration of the follower satellite is described by the gradient of
a potential function U such that
. U a (2.52)
Now we define an international terrestrial reference frame (ITRF). The ITRF, which is
fixed to the Earth and rotates with it, is used to describe the nonuniform gravity
perturbation depending on the position (the latitude and the longitude) of a satellite. Its
origin is at the center of the Earth, its first axis (
ˆ
ITRF
I ) extends through the point of
latitude 0º and longitude 0º, the second axis (
ˆ
ITRF
J ) is 90
to the east in the equatorial
plane, and the third axis (
ˆ
ITRF
K ) extends through the North Pole (Vallado [40]). Then, we
represent the potential function in the ITRF:
2
, , ,
2 1
,
1 sin
sin cos sin
.
l
l l ITRF
l ITRF
ITRF
l
l
ITRF
l m ITRF l m ITRF l m ITRF
l m ITRF
p ITRF
ITRF
R
J P
r
GM
U
r
R
P C m S m
r
GM
U
r
(2.53)
Here
l
J are the zonal harmonics coefficients, and
, l m
C and
, l m
S are the coefficients
associated with the tesseral and sectorial harmonics of the Earth, respectively; R
is the
equatorial radius of the Earth;
2 2 2
ITRF ITRF ITRF ITRF
r X Y Z is the distance from the
mass center of the Earth (which is generally assumed to coincide with the center of the
Earth) to the mass center of a satellite (where
ITRF ITRF ITRF ITRF
X Y Z
T
r is the
position vector of a satellite in the ITRF);
l
P is the Legendre polynomial of degree l ;
35
, l m
P is the associated Legendre function of degree l and order m ; and,
ITRF
and
ITRF
are the satellite’s latitude and longitude in the ITRF, respectively.
Next, let us consider the relative motion of satellites in the Hill frame. To begin with
we recall the two-body motion, that is,
2
.
body
GM
U
r
(2.54)
The corresponding acceleration represented in the Hill frame is given in Eq. (2.40) as
2 , 2
3/2
2
2 2
0 2 .
0
body Hill body
L L L L
L
U
r x r x r x r
GM
y y y
x r y z
z z z
a
RS RS
(2.55)
Here, the subscript “Hill” denotes that Eq. (2.55) is described in the Hill frame,
x y z
T
is the position vector of he center of mass of the follower satellite in the Hill
frame,
L
r is the distance from the center of the Earth to the leader satellite, and R is an
orthogonal rotation matrix that maps the ECI frame to the Hill frame, and each element of
the matrix R is given in Appendix A. The matrix S in Eq. (2.55) is the active rotation
matrix, which is the transpose of R .
Next, we consider the perturbing acceleration caused by the nonspherical shape of the
Earth. Substituting Eq. (2.53) (except for the first term
ITRF
GM
r
pertinent to two-body
attraction) into Eq. (2.52), we have (Vallado [40])
36
, , ,
, 2 2
2 2 2
1 1
,
p ITRF p ITRF p ITRF
ITRF
X ITRF ITRF ITRF
ITRF ITRF ITRF ITRF ITRF ITRF
ITRF ITRF ITRF
U U U
Z
d X Y
r r X Y
r X Y
(2.56a)
, , ,
, 2 2
2 2 2
1 1
,
p ITRF p ITRF p ITRF
ITRF
Y ITRF ITRF ITRF
ITRF ITRF ITRF ITRF ITRF ITRF
ITRF ITRF ITRF
U U U
Z
d Y X
r r X Y
r X Y
(2.56b)
2 2
, ,
, 2
1
,
p ITRF p ITRF ITRF ITRF
Z ITRF ITRF
ITRF ITRF ITRF ITRF
U U X Y
d Z
r r r
(2.56c)
where
, p ITRF
U is defined in Eq. (2.53), the position vector is given by
ITRF ITRF ITRF ITRF
X Y Z
T
r ,
ITRF ITRF
r r , and Eqs. (2.56) are evaluated in the ITRF.
Let C be the rotation matrix that maps the ITRF to the ECI frame, then we have
cos sin 0
sin cos 0 ,
0 0 1
t t
t t C
(2.57)
where
5
7.292115 10 rad/s
is the mean rotation rate of the Earth. Since R is the
rotation matrix that maps the ECI frame to the Hill frame, the resulting perturbing
acceleration described in the Hill frame is given by
,
, ,
,
X ITRF
aspherical Hill Y ITRF
Z ITRF
d
d
d
a RC
(2.58)
so that the acceleration of the follower expressed in the Hill frame is given by
2 , ,
,
, 3/2
2
2 2
,
:
= 0 2 ,
0
Hill body Hill aspherical Hill
L L L L X ITRF
Y ITRF
L
Z ITRF
r x r x r x r d
GM
y y y d
x r y z
z z z d
a a a
RS RS RC
(2.59)
37
where
, X ITRF
d ,
, Y ITRF
d , and
, Z ITRF
d are given in Eqs. (2.56). Here, it is noted that in Eq.
(2.59) the second derivative of the matrix S (or the matrix R) should be calculated, which
means that we need the time derivative of Eqs. (2.56) for the leader satellite. This
calculation is very complicated if done by hand, so the Maple code explicitly generating
this derivative is provided in Appendix B. In the code we can arbitrarily choose the
degree l .
2.2.5 Formation-keeping with attitude constraints
Finally, we include attitude constraints in formation-keeping. For example, it is important
to point solar panels of a satellite towards the Sun during the maneuver for stable power
supply. For attitude motion, quaternions are used to realize arbitrary orientations and
avoid singularities. Also, let us assume that the follower satellite has a diagonal inertia
matrix, J , in which the moments of inertia along its body-fixed principal axes of inertia
are placed. The body-fixed frame is used to describe the attitude motion of the follower
satellite. Each axis of the body-fixed frame of the follower satellite points along its
principal axes of inertia, and the origin of this frame is located at its center of mass. Here
we shall follow the method given in Udwadia and Schutte [42].
We begin with Lagrange’s equation
, , ,
u
d T T
t
dt
Γ u u
u u
(2.60)
38
where
0 1 2 3
u u u u
T
u is the quaternion 4-vector of the follower, , ,
u
t Γ u u is the
given generalized force vector, and T is the rotational kinetic energy of the follower
satellite which is given by
1
ˆ
.
2
T
T
ω J ω (2.61)
Here, the 4 by 4 augmented inertia matrix,
ˆ
J , is defined as
0
0
0 0 0
0 0 0
ˆ
: ,
0 0 0
0 0 0
x
y
z
J
J
J
J
J
0
J
0 J
(2.62)
where
0
J is an arbitrary positive number, and
x
J ,
y
J ,
z
J are the moments of inertia
along the principal axes of the follower satellite. Also, the 4 by 1 augmented angular
velocity vector, ω , is related with quaternions by
2 , 1,2, , , i N ω Eu (2.63)
where 0
x y z
T
ω
and the last three elements, described in the body frame,
are the angular velocities about the ECI frame of reference, and the 4 by 4 quaternion
matrix E is defined by
0 1 2 3
1 0 3 2
2 3 0 1
3 2 1 0
: .
u u u u
u u u u
u u u u
u u u u
T
1
u
E
E
(2.64)
Substituting Eq. (2.63) into Eq. (2.61) yields the kinetic energy in terms of quaternions
ˆ
2 . T
T T
u E JEu (2.65)
39
Then, assuming that there is no applied torque (i.e.
u
Γ 0 ), we now apply Lagrange’s
equation under the assumption that the components of the quaternion 4-vector are all
independent of each other, to obtain
ˆ ˆ ˆ
4 8 4 .
T T T
E JEu E JEu E JEu
(2.66)
This relation can be written in the form
u u
M u Q
by setting
ˆ
4
u
T
M E JE , and
u
Q to
be the right hand side of Eq. (2.66). It is noted that the mass matrix
ˆ
4
u
T
M E JE is
symmetric and positive definite, so it has always its inverse.
It is important to stress that up to now we have assumed that each component of the
quaternion vector u is independent of the others. However, to represent a physical
rotation of a rigid body we require the quaternion u to have unit norm so that
2
2 2 2 2
0 1 2 3
1. u u u u u (2.67)
After differentiating twice, we have the following constraint equation of the form of Eq.
(2.6)
0
1
0 1 2 3 0 1 2 3
2
3
u
u
u u u u u u u u
u
u
(2.68)
so that
0 1 2 3 0 1 2 3
, : .
u u
u u u u b u u u u N A u (2.69)
The resulting rotational equation of motion for the follower is thus given by the
fundamental equation, Eq. (2.10)
40
1
,
2
N
T -1
1
u E J ω Jω u u (2.70)
where
1
E , J , and N u are defined in Eqs. (2.64), (2.62), and (2.69), respectively, and
ω is a skew-symmetric matrix defined by
0
: 0 .
0
z y
z x
y x
ω
(2.71)
In Eq. (2.70), we have used the 8th-order system of differential equations using the
quaternion and its derivative 8-vector , u u instead of a 7th-order system that consists of
the quaternion 4-vector and the body-fixed angular velocity 3-vector , u ω . We can
readily switch from one to the other using Eq. (2.63) and its derivative, but we employed
in the current dissertation the straightforward Lagrangian formulation because of its
simplicity.
Finally, we combine the orbital and attitude dynamics. Defining the 7 by 1
generalized displacement vector t q as
0 1 2 3
, t x y z u u u u
T
q (2.72)
we have the following equation of unconstrained motion for the follower from Eqs.
(2.59) and (2.70)
, 1
2
Hill
t
N
-1
T -1
1
a M Q
a
E J ω Jω u u
(2.73)
41
where
Hill
a is a 3 by 1 vector on the right hand side of Eq. (2.59), and
1
2
N
T -1
1
E J ω Jω u u is a 4 by 1 vector on the right hand side of Eq. (2.70). Eq.
(2.73) describes the orbital and attitude dynamics of the follower satellite when no
constraint (control) is applied to the system.
When constraints (trajectory requirements) on either orbital and/or attitude motion,
which are of the form of Eq. (2.6), are applied, the generalized control force required to
follow the given constrained trajectory is explicitly determined by Eq. (2.8). In addition,
we can relate the 4 by 1 generalized quaternion torque,
u
Γ , which is determined by Eq.
(2.8), to the 3 by 1 physically applied torque about the body axis,
x y z
T
Γ
,
through the relation (Udwadia and Schutte [42])
0
1
.
2
u
EΓ
Γ
(2.74)
2.3 Numerical Example
In this section, an example is introduced to demonstrate the applicability of the control
methodology suggested in the previous sections. It is straightforward to extend this
example for application to more general situations. The numerical integration throughout
this paper is done in the Matlab environment, using a variable time step integrator with a
relative error tolerance of
8
10
and an absolute error tolerance of
12
10
.
42
We consider a system in which there are two follower satellites whose masses are
(1)
100 kg m and
(2)
50 kg m , where the superscripts “1” and “2” denote the first and
second follower satellites, respectively. Also, their moments of inertia along their
respective body-fixed axes are taken to be
(1) 2
400 200 150 kg m diag J and
(2) 2
200 100 75 kg m diag J , respectively. The leader satellite is in an orbit (in the
ECI frame), defined by
ˆ
ˆ ˆ cos , sin , sin ,
4
L
L L L L L L L L
n
X a n t Y a n t Z k t
(2.75)
where ˆ cos sin
4 4
L
L
n
a h R t
,
ˆ
cos
4
L
k R
,
5
6.2186 10 m h ,
6
6.3781 10 m R
is the equatorial radius of the Earth, and
3
1.0781 10 rad/s
L
n
.
Figure 2-2 shows the leader satellite’s somewhat complex trajectory in the ECI frame.
The trajectory in 3-dimensional space is given in the upper left figure. The upper right,
the lower left, and the lower right figures show the trajectories projected onto the XY ,
XZ , and YZ planes, respectively. In the lower left figure that depicts the projected
trajectory onto the XZ plane, the trajectory in the rear is invisible because it coincides
with the trajectory in the front. This orbit is chosen so as to demonstrate that the approach
developed herein is applicable even when the leader satellite’s orbit is arbitrary.
43
Fig. 2-2 The leader satellite’s trajectory in the ECI frame
The trajectory of the leader satellite is purposely chosen to be fairly complex. Its
orbital period is
4
2
4 2.3311 10 s 6.4754 hrs.
L
L
P
n
(2.76)
Its period in the XY plane is
4
L
P
and it moves between the latitudes of approximately
41
in the ECI frame. Additionally, its nominal altitude is a purposely taken to be a
44
periodic time-varying function which is approximately given by (we ignore the Earth’s
oblateness in the second member on the right hand side below)
2 2 2
2 2 2 2 2
1
cos sin 2 cos sin sin .
4 4 4 4 2 4
L
L L L
L L L
Altitude r R
X Y Z R
n n n
h R t hR t R t R
(2.77)
Figure 2-3 shows the somewhat complex desired time history of the altitude of the leader
satellite as it orbits around the Earth (solid line), and the corresponding desired variation
in its latitude (dotted line). We assume, without loss of generality, that time t is
measured from the instant that the ITRF coincides with the ECI frame, and we choose
2
L
P (twice the orbital period of the leader satellite) as the duration of time used for
numerical integration.
Fig. 2-3 Time history of the altitude and the latitude of the leader satellite
45
For the follower satellites we use the following orbital and attitude requirements: (1)
each follower satellite's orbit oscillates around a circle with radius when projected onto
the yz plane of the Hill frame with the leader satellite located at the center of the line
joining the two follower satellites; (2) each follower satellite, more specifically, the x
axis of the body frame of each follower satellite, points to the center of the Earth at all
times. Besides these two requirements, we shall impose the additional constraint that the
quaternion 4-vector of each follower satellite should have unit norm. These constraints
are summarized as
(1) (1) (1) (1)
2 , cos , sin , x z y t t z t t (2.78a)
(1)
(1) (1)
(1)
,
b
X
Y
Z
x T 0
(2.78b)
2 2 2 2
(1) (1) (1) (1) (1)
0 1 2 3
1, N u u u u u (2.78c)
(2) (2) (2) (2)
2 , cos , sin , x z y t t z t t (2.78d)
(2)
(2) (2)
(2)
,
b
X
Y
Z
x T 0
(2.78e)
and
2 2 2 2
(2) (2) (2) (2) (2)
0 1 2 3
1, N u u u u u (2.78f)
where
( ) i
x ,
( ) i
y ,
( ) i
z denote the position in the Hill frame of the i
th
follower satellite,
( ) i
X ,
( ) i
Y ,
( ) i
Z the position in the ECI frame of the i
th
follower satellite, the oscillation
46
amplitude t in Eqs. (2.78a) and (2.78d) is chosen as
0
sin 0.1 0.025 sin 4
L L L
t t a a n t where
6
7.0 10 m
L
a is the
maximum value of ˆ cos sin
4 4
L
L
n
a h R t
or the value of ˆ
L
a at the initial time, i.e.
each follower satellite oscillates sixteen times (in the Hill frame) about its nominal
circular orbit (of radius
0
0.1
L
a ) during one period (
L
P ) of the leader satellite’s
motion. The rotational frequency in Eqs. (2.78a) and (2.78d) is set to equal
L
n , and the
follower satellites, in addition, oscillate four times as they move around their nominal
circular orbits (as seen in the Hill Frame) for each revolution they make around the leader
satellite. The position vector in the ECI frame in Eqs. (2.78b) and (2.78e) can be
transformed to the one in the Hill frame, and vice versa, by using the relation
( ) ( )
( ) ( )
( ) ( )
,
i i
L
i i
i i
X x r
Y y
Z z
T
R
(2.79)
where the components of the transformation matrix R are given in Appendix A. In Eqs.
(2.78b) and (2.78e),
b
x is the skew-symmetric matrix given by
0 0 0
0 0 1
0 1 0
b
x
(2.80)
corresponding to the unit vector along the x axis of the body frame 1 0 0
b
T
x ,
and
( ) i
T is a transformation matrix that maps the ECI frame into the body frame of the i
th
follower satellite which is of the form
47
2 2 2 2
2 2 2 2
2
( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )
0 1 2 3 1 2 0 3 0 2 1 3
( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )
1 2 0 3 0 1 2 3 2 3 0 1
( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )
1 3 0 2 0 1 2 3 0 1
2 2
2 2
2 2
i i i i i i i i i i i i
i i i i i i i i i i i i i
i i i i i i i i i
u u u u u u u u u u u u
u u u u u u u u u u u u
u u u u u u u u u u
T
2 2 2
( ) ( ) ( )
2 3
.
i i i
u u
(2.81)
Eqs. (2.78b) and (2.78e) originate from the fact that the desired pointing axis (i.e. x
axis of the body frame) is constrained to point along the vector connecting each follower
satellite and the center of the Earth,
( ) ( ) ( ) i i i
X Y Z
T
. The components of this
vector, in turn, are transformed into the body frame by the transformation matrix
( ) i
T ,
and the cross product of this transformed vector and the x axis of the body frame is
zero because they are parallel.
We choose the initial conditions for orbital motion of the two follower satellites as
(1) 4 (1) 5 (1) 4
0 1.7493 10 m, 0 6.9300 10 m, 0 3.4111 10 m, x y z
(1) (1) (1)
0 377.3473 m/s, 0 0 m/s, 0 754.6946 m/s, x y z (2.82a)
(2) 4 (2) 5 (2) 4
0 1.7493 10 m, 0 6.9300 10 m, 0 3.4111 10 m, x y z
(2) (2) (2)
0 377.3473 m/s, 0 0 m/s, 0 754.6946 m/s, x y z (2.82b)
and the initial conditions for attitude motion as
(1) (1) (1) (1)
0 1 2 3
0 0.001924733, 0 0.002997640, 0 0.2, 0 0.9797894, u u u u
(1) (1) (1) (1)
0 1 2 3
0 0, 0 0, 0 0, 0 0, u u u u (2.83a)
(2) (2) (2) (2)
0 1 2 3
0 0.001924733, 0 0.002997640, 0 0.2, 0 0.9797894, u u u u
(2) (2) (2) (2)
0 1 2 3
0 0, 0 0, 0 0, 0 0. u u u u (2.83b)
48
The two follower satellites have the same initial components of their quaternions. It must
be noted that the initial conditions given by Eqs. (2.82) and (2.83) do not satisfy the
constraint Eqs. (2.78) (except Eqs. (2.78c) and (2.78f), since otherwise we would not be
describing a physical rotation), so we should employ the modified constraint equation,
Eq. (2.48). These initial discrepancies are a consequence of errors in the initial insertion
of the follower satellites into their desired formation configuration. The constraint
equation can be compactly written in vector form as
( ) ( ) ( )
, 1,2,
i i i
i Φ Φ Φ 0
(2.84)
where
( ) i
Φ is a 6 by 1 vector devised from the constraints, Eqs. (2.78)
2 2 2 2
( ) ( )
( )
( )
( )
( ) ( ) ( )
21 22 23
( ) ( ) ( )
31 32 33
( ) ( ) ( ) ( )
0 1 2 3
2
cos
sin
.
1
i i
i
i
i
i i i
i i i
i i i
i i i
i i i i
x z
y t t
z t t
T X T Y T Z
T X T Y T Z
u u u u
Φ
(2.85)
In this paper, 0.001 and 0.0005 are chosen for all the constraints, and by setting
( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )
0 1 2 3
i i i i i i i i
t x y z u u u u
T
q
, Eq. (2.84) is rewritten in the form of Eq.
(2.6) to obtain A and b which are used for determining the exact control force by Eq.
(2.8). Also, the fourth and fifth components of
( ) i
Φ can be rewritten with
( ) i
x ,
( ) i
y , and
( ) i
z in the Hill frame using Eq. (2.79).
49
Figures 2-4 and 2-5 represent the orbits of the first and second follower satellites with
control projected on the yz plane and xz plane in the Hill frame, respectively. The
scale is normalized by
0
0.1
L
a . As seen in the figures, in these normalized coordinates
each follower satellite oscillates around the circumference of a circle of radius unity
(shown by the dotted line) four times during each revolution around the leader satellite.
The figures also show that the controlled motion of the follower satellites asymptotically
approach the desired trajectories (described by Eqs. (2.78a) and (2.78d)) though they start
off with initial orbit insertion errors. Figures 6 and 7 represent the orbits of the first and
second follower satellites without control projected on the yz plane and xz plane in
the Hill frame, respectively. As seen from the figures, without control, the follower
satellites show complex motions, nothing like those in Figs. 2-4 and 2-5. In Figs. 2-8 and
2-9, the time history of each component of the quaternion for both follower satellites is
shown where time is normalized by
L
P , the period of the leader satellite (see Eq. (2.76)).
Because they have the same initial components of quaternion (see Eqs. (2.83)), the two
figures (Figs. 2-8 and 2-9) are similar at the beginning but differ from one another as time
goes by.
50
Fig. 2-4 Controlled motion in the yz- and xz-planes of the first follower satellite. The follower
satellite oscillates 4 times about the circle of radius
0
for each revolution around the leader satellite
(as seen in the Hill Frame)
Fig. 2-5 Controlled motion in the yz- and xz-planes of the second follower satellite
51
Fig. 2-6 Uncontrolled motion in the yz- and xz-planes of the first follower satellite
Fig. 2-7 Uncontrolled motion in the yz- and xz-planes of the second follower satellite
52
Fig. 2-8 Quaternion of the first follower satellite
Fig. 2-9 Quaternion of the second follower satellite
In Figs. 2-10 and 2-11 we show the explicitly obtained control forces per unit mass of
the first and second follower satellites to maintain the desired formation and their total
53
magnitude. The force components are described in the Hill frame, and time is normalized
to the period of the leader satellite (i.e.
L
P ). As seen, relatively large control forces are
brought into play initially in order to eliminate the large initial errors in the orbits of the
two follower satellites; also, the complex orbit of the leader satellite contributes to an
increased required thrust force. Although it seems impractical to execute such large-thrust
maneuvers initially, one can reduce magnitude of the control thrust and torque by
choosing smaller
i
’s and
i
’s than those used here; but this would, of course, come at
the expense of a longer time duration needed to get the follower satellites to their
required trajectories. Figures 2-12 and 2-13 illustrate the control torques per unit mass of
the first and second follower satellites for attitude requirements. The torque components
are described in the body frame of each follower satellite, and it is seen again that
relatively large control torques are required to mitigate initial attitude errors. This is
because although the initial quaternions satisfy the constraints (unit norm), their
attitudinal motion is dependent on their orbital motion, thus initial errors in orbital motion
yield errors in attitudinal motion. After these initial adjustments, as seen from the figures,
the torques required for attitude control become very small. Figures 2-14 and 2-15
represent errors in the satisfaction of the desired trajectory requirements described by Eq.
(2.85). We denote these errors by
( ) ( ) ( )
1
2
i i i
e t x z ,
( ) ( )
2
cos
i i
e t y t t ,
( ) ( )
3
sin
i i
e t z t t ,
( ) ( ) ( ) ( )
4 21 22 23
i i i i
i i i
e t T X T Y T Z ,
54
( ) ( ) ( ) ( )
5 31 32 33
i i i i
i i i
e t T X T Y T Z , and
2 2 2 2
( ) ( ) ( ) ( ) ( )
6 0 1 2 3
1
i i i i i
e t u u u u , where the
superscripts 1,2 i refer to the two follower satellites. We see that initially these errors
are very large (except for
( )
6
i
e ) since the follower satellites have not been inserted into
orbit with the correct initial conditions that satisfy the desired trajectory requirements,
however, they converge to zero as time progresses. More specifically, for the first
follower satellite, these initial errors
(1)
0
j
e , 1,2,3,4,5,6 j are of the order of
2
10 meter ,
3
10 meter ,
4
10 meter ,
5
10 meter ,
5
10 meter , and
12
10
, respectively.
However, they converge to zero as time progresses, and at the end of the duration of
integration, 2
L
P (twice the period of the leader satellite), the final errors
(1)
2
j L
e P ,
1,2,3,4,5,6 j are of the order of
8
10 meter
,
5
10 meter
,
5
10 meter
,
5
10 meter
,
5
10 meter
, and
12
10
, respectively. Figures 2-16 and 2-17 verify this. As seen from the
time histories over the duration 1.95 2
L L
P P , the six errors become very small compared
with the initial values (see Figs. 2-14 and 2-15).
55
Fig. 2-10 Required control forces per unit mass in the Hill frame and the force magnitude for
the first follower satellite
Fig. 2-11 Required control forces per unit mass in the Hill frame and the force magnitude for
the second follower satellite
56
Fig. 2-12 Required control torques per unit mass in the body frame and the torque magnitude
for the first follower satellite
Fig. 2-13 Required control torques per unit mass in the body frame and the torque magnitude
for the second follower satellite
57
Fig. 2-14 Errors for the first follower satellite in the satisfaction of the six trajectory
requirements
Fig. 2-15 Errors for the second follower satellite in the satisfaction of the six trajectory
requirements
58
Fig. 2-16 Errors in the first follower satellite’s trajectory requirements over the time duration
1.95 2
L L
P P
Fig. 2-17 Errors in the second follower satellite’s trajectory requirements over the time duration
1.95 2
L L
P P
59
Chapter 3
SEPARATION PRINCIPLE: EXPLICIT
DETERMINATION OF THE CONTROL FORCE AND
TORQUE FOR SATELLITE FORMATION-KEEPING
3.1 Determining the Explicit Form of the Control Forces and Torques
We begin this chapter by recalling the following equation of unconstrained motion for
each follower satellite in the ECI frame:
1
( ) ( )
( )
2
( ) ( ) ( )
( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )
: ,
1
2
i i
i
body aspherical
i i i
i i i i i i i i
t
N
T
-1
1
2 1
a a
a
a M Q
a E J ω J ω u u
(3.1)
where the superscript “i” is pertinent to a quantity of the ith follower satellite, and the
upper block in Eq. (3.1) is a 3-vector related to the orbital part:
2 2 2
( ) ( )
,
( ) ( ) ( ) ( ) ( )
2 , 3/2
( ) ( ) ( )
( ) ( )
,
,
i i
X ITRF
i i i i i
body aspherical Y ITRF
i i i
i i
Z ITRF
X d
GM
Y d
X Y Z
Z d
1
a a a C
(3.2)
where
( ) ( ) ( ) i i i
X Y Z
T
is the position vector of the ith follower satellite in the ECI
frame, and the matrix C and
( ) ( ) ( )
, , ,
i i i
X ITRF Y ITRF Z ITRF
d d d
T
are given in Eqs. (2.57) and
(2.56), respectively. Also, the lower block in Eq. (3.1) is a 4-vector related to the attitude
part:
1
( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )
1
,
2
i i i i i i i i
N
T
2 1
a E J ω J ω u u
(3.3)
60
which is given in Eq. (2.70). Also, the 7 by 7 positive definite mass matrix of the ith
follower satellite
( ) i
M is defined as
( )
( )
( )
: ,
i
i
i
m
3 3 3 4
4 3 u
I 0
M
0 M
(3.4)
where
( ) i
m is the mass of the ith follower satellite,
3 3
I is the 3 by 3 identity matrix, and
( ) ( ) ( ) ( )
ˆ
4
i i i i
T
u
M E J E is a positive definite 4 by 4 matrix, where
( ) i
E and
( )
ˆ
i
J are defined in
Eqs. (2.64) and (2.62). Then Eq. (3.1) describes the orbital and attitude dynamics of the
ith follower satellite when no constraint (control) is applied to the system.
Let
( ) i
t
C
F be the control force 3-vector for orbital control and
( ) i
t
C
u
Γ be the
generalized quaternion torque 4-vector for attitude control, then they form the constraint
(control) force 7-vector
( ) i
t
C
Q in Eqs. (2.8), that is,
( )
( )
( )
.
i
i
i
t
t
t
C
C
C
u
F
Q
Γ
(3.5)
When constraints on either orbital and/or attitude motion, which are of the form of Eq.
(2.6), are applied, the constraint force and torque
( ) i
t
C
Q required to follow the given
constrained trajectory are explicitly determined by Eqs. (2.8). In addition, one could
relate the generalized quaternion torque 4-vector,
( ) i C
u
Γ to the physically applied torque 3-
vector about the body axis,
( ) ( ) ( ) ( ) i i i i
x y z
T
C
Γ
, through the relation Eq. (2.74) or
( ) ( )
1
, 1,2, , ,
2
i i
i N
C C
1 u
Γ E Γ (3.6)
where
1
E is defined in Eq. (2.64).
61
With the aid of the analytical expression for the control forces (
( ) i
t
C
F ) and torques
(
( ) i
t
C
u
Γ ) given by Eqs. (2.8), it is possible to obtain these control functions in closed
form once orbital and/or attitude constraints are given. These analytical solutions provide
many important physical insights into the dynamics of multiple satellites, which will be
shown in the next section that cannot be easily inferred from numerical computations.
The methodology adopted in this chapter is general in the sense that the explicit form of
the control forces and torques is immediately obtained once constraints are given. Now,
the following three kinds of constraints are considered in this chapter.
Let us first assume that for the ith follower satellite the following form of the
constraints is obtained from s consistent orbital constraints:
, , , , , t t
11 1
A X X X b X X
(3.7)
where
11
A is an s by 3 matrix, X Y Z
T
X
is the acceleration 3-vector of the ith
follower satellite in the ECI frame,
1
b is an s-vector, and s 3 is the number of the
given (independent) orbital constraints. Also, as in Eq. (3.7) the superscript “i” will be
omitted for brevity from here on, unless required for clarity.
As a specific example, let us consider the projected-circular orbit (PCO) exemplified
in Eqs. (2.12a) and (2.12b):
2 2 2
, y z (3.8a)
2 , x z (3.8b)
62
where there are two constraints, so s is 2. Eq. (3.8a) states that the follower satellite
should move on a circle with constant radius , where is the constant distance
between the leader satellite and the ith follower satellite when projected on the local
horizontal y z plane corresponding to it. Eq. (3.8b) makes the relative motion
bounded for every axis, and it also matches the solutions of the HCW equations
satisfying the constraints Eq. (3.8a) (Sabol et al. [43]). Each follower satellite can have its
own orbital requirements with respect to the leader satellite that can be correspondingly
expressed as, for example, in Eqs. (3.8). By transforming Eqs. (3.8) to the ECI frame the
two orbital constraints can be expressed as,
2 2
2
,
2 3
R X R X (3.9a)
2 ,
L
r
1 3
R X R X (3.9b)
where
1
R ,
2
R , and
3
R are the first, second, and third row vectors of the rotation matrix
R in Eq. (2.30), respectively, and they are explicitly given in Appendix A. To obtain a
general form of the constraint equation expressed by Eq. (2.6), the following are defined:
2 2
2
1
: ,
2 3
R X R X (3.10a)
2
: 2 ,
L
r
1 3
R X R X (3.10b)
so that
1
2 2 ,
2 2 2 3 3 3
R X R X R X R X R X R X
(3.11a)
2
2 ,
L
r
1 1 3 3
R X R X R X R X
(3.11b)
63
and
2
1
2
2 2 2
2 2 2 ,
2 2 2 2 2 2
3 3 3 3 3 3
R X R X R X R X R X R X
R X R X R X R X R X R X
(3.12a)
2
2 2 2 .
L
r
1 1 1 3 3 3
R X R X R X R X R X R X
(3.12b)
Then, using Eq. (2.48) one can obtain the constraint equations in the form of Eq.
(2.6):
2 2
1 1 1 1
2 2 2 2 2
2 2 ,
2 2 3 3 2 2 2 2 2 3 3
3 3 3
R X R R X R X R X R X R X R X R X R X R X
R X R X R X
(3.13a)
2 2 2 2
2 2 2 2 2 .
L
r
1 3 1 1 3 3
R R X R X R X R X R X
(3.13b)
Eqs. (3.13) can be compactly written as
2 2
1 1 1 1
2 2 2 2
2 2 2 2
2
2 2
,
2
2 2 2 2
L
r
2 2 2 2 2 3 3
2 2 3 3
3 3 3
1 3
1 1 3 3
R X R X R X R X R X R X R X
R X R R X R
R X R X R X
X
R R
R X R X R X R X
(3.14)
which is of the form of Eq. (3.7) with
2
, , ,
2
t
2 2 3 3
11
1 3
R X R R X R
A X X
R R
(3.15a)
2 2
1 1 1 1
2 2 2 2
2 2 2 2
2 2
, , ,
2 2 2 2
L
t
r
2 2 2 2 2 3 3
3 3 3
1
1 1 3 3
R X R X R X R X R X R X R X
R X R X R X
b X X
R X R X R X R X
(3.15b)
64
where
11
A is a 2 by 3 matrix and
1
b is a 2-vector. It is noted in passing that the matrix
11
A and the vector
1
b do not involve any quaternions. As with Eqs. (3.8a) and (3.8b), it
shall be assumed from here on that the 3 s orbital constraints for each of the follower
satellites are independent of one another, so that the rank of the matrix
11
A is full.
Next, let us consider attitude constraints for target tracking. This constraint has been
already handled in the previous chapter, but here it shall be more rigorously investigated
in this chapter. It is assumed that the z-axis of the body frame of the ith follower satellite
described by the unit vector
( )
ˆ 0 0 1
i
T
body
z is required to point to a given target and
that the target’s position vector in the ECI frame, which is given by
( ) ( ) ( ) ( ) i i i i
t t t
t X t Y t Z t
T
t
X , may be an arbitrarily prescribed time-varying
function. This constraint can be mathematically expressed using the fact that in the body
frame of reference the vector
( )
ˆ
i
body
z should be parallel to the vector
( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )
:
i i i i i i i i i
t t t
t t t X t X t Y t Y t Z t Z t
T
t
ξ X X 0
connecting the follower satellite and the target. As before, the superscript “i” shall be
dropped in what follows, for brevity, but note that each follower satellite is required to
track, in general, a different target whose position may be arbitrarily prescribed in time.
The aim is to perform, along with precise orbital control, attitude control so that each
follower satellite can precisely track its assigned target.
65
Let us recall the transformation matrix T that maps the ECI frame into the body
frame of the ith follower satellite which is given in Eq. (2.81):
2 2 2 2
11 12 13 0 1 2 3 1 2 0 3 1 3 0 2
2 2 2 2
21 22 23 1 2 0 3 0 1 2 3 2 3 0 1
2 2 2 2
31 32 33 1 3 0 2 2 3 0 1 0 1 2 3
2 2
: : 2 2 ,
2 2
T T T u u u u u u u u u u u u
T T T u u u u u u u u u u u u
T T T u u u u u u u u u u u u
1
2
3
T
T T
T
(3.16)
where
1
T ,
2
T , and
3
T are the first, second, and third row vectors of the T matrix,
respectively. Then t Tξ is the vector connecting the follower satellite and the target
expressed along the body frame, and the attitude constraint for target pointing becomes
ˆ ˆ : .
1
body body 2
3
T ξ
z Tξ z T ξ 0
T ξ
(3.17)
Eq. (3.17) describes the attitude tracking requirement for the ith follower satellite, and
each of the N follower satellites can have their own time-varying attitude requirements
for tracking different targets.
Expanding Eq. (3.17) for the ith follower satellite, one has
11 12 13
21 22 23
31 32 33
0 1 0 0 1 0 0
1 0 0 1 0 0 0 .
0 0 0 0 0 0 0
t t t
t t t
t t t
T X X T Y Y T Z Z
T X X T Y Y T Z Z
T X X T Y Y T Z Z
1
2
3
T ξ
T ξ
T ξ
(3.18)
Then, the two constraints for target pointing are given by
21 22 23
0,
t t t
T X X T Y Y T Z Z
(3.19a)
11 12 13
0.
t t t
T X X T Y Y T Z Z (3.19b)
Eqs. (3.19a) and (3.19b) can also be compactly written, for future use, as
66
0,
2
T ξ (3.19c)
0.
1
T ξ (3.19d)
As before, the following are defined:
3 3
: , ,
2 2 2
T ξ T ξ T ξ
(3.20a)
4 4
: , ,
1 1 1
T ξ T ξ T ξ
(3.20b)
so that
3 4
0, 0. (3.21)
Remark: Instead of Eq. (3.17), one might think that we can use only one constraint
equation like the following:
ˆ ˆ ˆ cos ,
T T
body body body
z Tξ z Tξ z Tξ ξ (3.22)
because ˆ
body
z is a unit vector, T is an orthogonal matrix, and , which is the angle
between the two vectors ˆ
body
z and T , is zero. Eq. (3.22) can be rewritten as
1/2
ˆ ,
T T
body 3
z Tξ T ξ ξ ξ ξ (3.23)
where
3
T is the third row vector of the matrix T . Because T is an orthogonal matrix,
from Eq. (3.23) we have
1/2 1/2 1/2 1/2
2 2 2
T T T T T T
3 3 3 1 1 2 2 3 3 1 2 3
T ξ ξ ξ ξ I ξ ξ T T T T T T ξ T ξ T ξ T ξ (3.24)
or simply
1/2
2 2 2
,
3 1 2 3
T ξ Tξ T ξ T ξ (3.25)
67
which is equivalent to saying that
2 2
2 2
3 4
0.
1 2
T ξ T ξ (3.26)
Eq. (3.26) represents a single constraint rather than the existing two constraints Eq.
(3.21). Mathematically Eq. (3.26) and (3.21) are equivalent if
3
and
4
have real values.
However, the single constraint Eq. (3.26) cannot be used as the target-pointing attitude
constraint due to the following reason. Let us assume that differentiating Eq. (3.21) twice
with respect to time are
3
0
and
4
0
, and they yield
, ,
3 3 4 4
A q b A q b
(3.27)
whereas differentiating Eq. (3.26) twice with respect to time yields
2 2
3 3 4 4 3 4
. (3.28)
However, since
3
0 and
4
0 , Eq. (3.28) just results in
0
1 7
0 q
(3.29)
and Eq. (3.29) cannot be used as the constraint equation because any q satisfies Eq.
(3.29). In brief, Eq. (3.26) should not be used as the target-pointing attitude constraint
because it does not provide any information about constraints.
Besides these two requirements Eqs. (3.19c) and (3.19d), it is necessary to impose an
additional constraint, as before, that the quaternion 4-vector of each follower satellite
68
must have unit norm, so that physical rotations are described. This constraint is of the
form
2 2 2 2
0 1 2 3
: 1, N u u u u u (3.30)
so that
2 2 2 2
5 0 1 2 3 5 0 0 1 1 2 2 3 3
: 1, 2 . u u u u u u u u u u u u (3.31)
Following the same procedure in the orbital part, one can have the following attitude
constraint equations:
,
21 22 2
X
A A b
u
(3.32)
where
21 22 23
11 12 13
: ,
0 0 0
T T T
T T T
21
21
1 3
A u
A u
0
(3.33)
3 0 1 2 1 0 1 2 3 0 3 2
0 3 2 1 2 3 2 1
2 2 2 2
, , 2 2 2
t t t t
t t t t
t t t t
t t
t t
t t
X X X X X X X X
u u u Y Y u u u Y Y u u u Y Y u u u Y Y
Z Z Z Z Z Z Z Z
X X X X
u u u Y Y u u u Y Y u u
Z Z Z Z
22 t
A X X u
0 3 0 1
0 1 2 3
2
2 2 2 2
, ,
:= ,
2
t t
t t
t t
X X X X
u Y Y u u u Y Y
Z Z Z Z
u u u u
22 t
T
A X X u
u
(3.34)
69
2 2 2 2
1 2 0 3 0 1 2 3 2 3 0 1
3 3 3 3
21 22 23 21 22 23
2 2 2 2
0 1 2 3 1 2
4 2 4
2 2 2
2 4
, , , , , ,
t t t
t t t t t t
t
u u u u X X u u u u Y Y u u u u Z Z
T X X T Y Y T Z Z T X T Y T Z
u u u u X X u u
2 t t t
b X X u X X u X
0 3 1 3 0 2
4 4 4 4
11 12 13 11 12 13
5 5 5 5
4
2 2 2
2
t t
t t t t t t
u u Y Y u u u u Z Z
T X X T Y Y T Z Z T X T Y T Z
N
u
5 5 5 5
, , , , , ,
: ,
2N
2 t t t
b X X u X X u X
u
(3.35)
where
2 2 2 2
0 1 2 3
: N u u u u u . The matrix
21
A is a 3 by 3 matrix,
22
A is a 3 by 4
matrix, and
2
b is a 3-vector; to obtain the corresponding barred quantities, the last row in
each of these matrices needs to be dropped. In Appendix C it is shown that
22
A always
has full rank.
Finally, the constraint equation which is of the form of Eq. (2.6), including both
orbital and attitude constraints, is written as:
: : : ,
11 1
21 22 2
A 0 b X
Aq b
A A b u
(3.36)
where the matrix A is in general (s+3) by 7 and 3 s is the number of independent
orbital constraints or the number of the independent rows of the matrix
11
A . Then, the
control force and torque are explicitly calculated using Eq. (2.8). In Appendix C, it is
proved that the matrix A in Eq. (3.36) always has full rank. Then, one can use the
regular inverse instead of the MP-generalized inverse of the matrix
-1 T
AM A in Eq. (2.8)
so that the control forces ( t
C
F ) and torques ( t
C
u
Γ ) on the ith follower satellite are
given explicitly by the 7-vector
70
: .
t
t
t
C
-1
C T -1 T
C
u
F
Q A AM A b Aa
(3.37)
Using the regular inverse as in Eq. (3.37) instead of the MP-generalized inverse as in Eq.
(2.8) saves considerable computation time and dramatically improves numerical
accuracy. The first component of t
C
Q which is a 3-vector, t
C
F , is the control force
for orbital control, whereas the second component which is a 4-vector, t
C
u
, is the
generalized control torque for attitude control. In Appendix D it is shown that the right
hand side of Eq. (3.37) simplifies to yield
, t m m m m
C T T
O 1 21 2 21 1 22 2 21 O O 1 21 O,A
F α P A D b A a A a A α α P A α (3.38)
, t m m
C T T
u 22 2 21 1 22 2 21 O 22 O,A
A D b A a A a A α A α (3.39)
where
, , : , , , : , , , t t t
-1 -1
T T + T T +
0 11 11 11 11 11 11 1 3 3 0 3 3 11 11 11 11 3 3 11 11
P X X A A A A A A P X X I P X X I A A A A I A A
(3.40)
, , : , , , , , : , t t m
-1 -1
T T + T -1 T
O 11 11 11 1 11 1 11 1 11 1 t 21 1 21 22 u 22
α X X A A A b A a A b A a D X X u X A P A A M A
(3.41)
and
, , , , , , , : . t
O,A t t t 2 21 1 22 2 21 O
α X X u X X u X D b A a A a A α
(3.42)
The matrices
0
P and
1
P are orthogonal projection operators because they are both
symmetric idempotent matrices. The matrix
0
P is a projection on to the column space of
T
11
A , and the matrix
1
P is the projection on to the null space of
11
A . It should be noted
that the matrix D in Eq. (3.41) is well defined since the 3 by 4 matrix
22
A has full rank
(see Appendix C) so that
-1 T
22 u 22
A M A is positive definite, and
T
21 1 21
A P A is semi-positive
71
definite since
1
P is idempotent. Let us briefly investigate the relationship between
0
P and
1
P . First it should be pointed out that (Udwadia and Kalaba [44])
. rank rank rank s
+
0 11 11 11
P A A A (3.43)
Since the column space of
+
1 3 3 11 11
P I A A is the same as the null space of
11
A
(Udwadia and Kalaba [44]), the rank of
1
P is
3 . rank nullity s
1 11
P A (3.44)
Combining Eqs. (3.43) and (3.44) yields
3 rank rank
0 1
P P (3.45)
and Eqs. (3.44) and (3.45) will be used later.
The structure of Eqs. (3.38) and (3.39) is better exposed when written in the form
,
1
2
t
m
t
T
C
3 3 1 21
O
C
T
O,A
1 22
I P A
α F
α Γ
0 E A
(3.46)
where Eq. (3.6) has been used in Eq. (3.46) to obtain the physical torque 3-vector, t
C
Γ ,
applied to the follower satellite. Eq. (3.46) is valid whether or not there is sensor
measurement noise that causes the initial conditions, both in orbit and attitude, of the
follower satellites to be incorrect at the start of the maneuver, resulting in initial insertion
errors.
One can further particularize the general result given in Eq. (3.46) when the attitude
constraints Eqs. (3.19) and the unit-norm constraint Eq. (3.30) are satisfied and there are
72
no attitude sensor measurement errors. As pointed out earlier, the attitude constraints will
be asymptotically satisfied in the presence of sensor measurement noise. Recalling that
(see Eqs. (3.33)-(3.35)) if the attitude constraints and the unit-norm constraint are met,
, , ,
2 2 N
2 21 22
T 21 22 2
1 3
b A A
A A b
u 0 u
(3.47)
and one can find that the control force and the physically applied torque can be rewritten
as (see Appendix E):
, t m m
C T
O 1 21 O,A
F α P A α (3.48)
,
2
m
t
C T
1 22 O,A
Γ E A α (3.49)
where
, , , , , , , , , , , , ,
4
m
t t
-1
T T -1 T
t 21 1 21 22 1 1 22 O,A t t t 2 21 1 22 2 21 O
D X X u X A P A A E J E A α X X u X X u X D b A a A a A α
(3.50)
with
1
P and
O
α defined in Eqs. (3.40) and (3.41). Also, the invertibility of the matrix D
in Eq. (3.50) is guaranteed. Note that the matrix
-1 T
22 u 22
A M A is positive definite where
0
0
0
1
1
2
4 2
1 1
2
4 2
1
0,
4
4
J
J
J
T
-1 T T T 22
22 u 22 1 22
T
-1 1
T T -1 T 22
1 1 22
T
T -1 T
22 1 1 22
0
u A
A M A u E A u
E u
0 J
A
uu E J E A u
u
A E J E A 0
0
(3.51)
and Eqs. (2.62), (2.64), (3.47), and the identity
1
E u 0 have been used in deriving Eq.
(3.51), assuming that the attitude constraints and the unit-norm constraint, Eqs. (3.19c),
73
(3.19d), and (3.30), are satisfied so that
22
A u 0 in Eq. (C-3) of Appendix C. From Eq.
(3.51) it is concluded that the matrix
4
m
T T -1 T
21 1 21 22 1 1 22
A P A A E J E A is positive definite
since
T -1 T
22 1 1 22
A E J E A is positive definite and
T
21 1 21
A P A is semi-positive definite so that D
(see Eq. (3.50)) is well defined. Eq. (3.46) now reads
,
,
1
2
t
m
t
T
C
3 3 1 21
O
C
T
O A
1 22
I P A
α F
α Γ
0 E A
(3.52)
where the barred quantities are given in Eqs. (3.33)-(3.35),
O,A
α is defined in Eq. (3.50),
and
1
P in Eq. (3.40).
In conclusion, Eqs. (3.46) and (3.52) that explicitly give the control forces and control
torques form part of the central results of this chapter. It is important to note that Eq.
(3.46) is more general than Eq. (3.52) in that Eq. (3.46) is applicable to cases when the
attitude constraints Eqs. (3.19c) and (3.19d) are not exactly satisfied, for example, when
initial condition errors caused by sensor measurement noise are present. Eq. (3.52) is
simpler than Eq. (3.46), but it is obtained based on the assumption that the attitude
constraints Eqs. (3.19c) and (3.19d) are satisfied throughout the maneuver (including at
the time of insertion), though not necessarily the orbital constraints. More precisely, Eq.
(3.52) requires the unit-norm constraint Eq. (3.30) as well. If this constraint is not
satisfied, it is impossible to describe physical rotation; hence the quaternions have unit
74
norm at the initial time and are controlled to maintain this norm by the controller
developed in this chapter.
Also, the acceleration vector of the constrained (controlled) system is explicitly
obtained from Eqs. (3.46) as
,
,
,
m
m
m
-1
T
O 3 3 1 3 3 1 21
T
O A u 2 22
T
O 1 3 3 1 21
-1 T
O A 2 u 22
α I 0 a I P A X
q
α 0 M a 0 A u
α a I P A
α a 0 M A
(3.53)
or from Eqs. (3.52) as
,
,
.
m
m
m
-1
T
O 3 3 1 3 3 1 21
T
O A u 2 22
T
O 1 3 3 1 21
-1 T
O A 2 u 22
α I 0 a I P A X
q
α 0 M a 0 A u
α a I P A
α a 0 M A
(3.54)
As before, Eq. (3.54) is restricted to situations when the attitude constraints Eqs. (3.19c)
and (3.19d) (and the unit-norm constraint Eq. (3.30)) are satisfied, Eq. (3.53) being the
more general result. From Eqs. (3.41) and (3.42) it is seen that
O
α is a function only of
the orbital constraints (requirements) while
, O A
α depends on both the orbital and the
attitude constraints (requirements). Thus from Eq. (3.46) when
T
1 21 O,A
P A α 0, the control
force t
C
F required to be applied to the follower satellite for orbital control does not
depend on the attitude of the follower satellite or the attitude constraints it is required to
satisfy. On the other hand, the control torques, t
C
Γ , are dependent, in general, on both
75
the orbital and the attitude requirements, and they depend on the orbital position and
velocity of the follower satellite, and of course on u and u . This leads to an important
principle explained in the next section.
3.2 Separation Principle
The structure of Eq. (3.46) leads to the following result.
Separation Principle. When the independent orbital trajectory constraints (trajectory
requirements) for the ith follower are such that
, , , t t
11 1
A X, X X b X, X
(3.55)
where t X t Y t Z t
T
X is the position vector of the ith follower satellite in
the ECI frame, and the attitude pointing constraints (requirements) are such that
,
21 22 2
X
A A b
u
(3.56)
where u is the quaternion 4-vector describing the orientation of the follower satellite,
then the control force t
C
F on the follower satellite needed at time t to satisfy the
orbital constraints (requirements) will be independent of any attitude pointing constraints
(requirements) at time t provided at that time
T
1 21 O,A
P A α 0 (3.57)
is satisfied, where
1
P and
O,A
α are respectively defined in Eqs. (3.40) and (3.42). The
corresponding control force at time t is given by , t m
C
O
F α where
O
α is explicitly
76
given in Eq. (3.41), and it does not depend on the attitude of the follower satellite. Eq.
(3.57) can be often satisfied in one of two different ways, leading to the following two
cases.
Case 1: When
1
P 0 .
For any controller that does exact orbital control, for all time, is independent of attitude
control (or any attitude consideration) if
11
A is invertible (i.e. if there are three
independent orbital constraints). Only these three orbital constraints fully determine the
orbital state vector t X and t X
with the given initial conditions. Regarding the
controller employed in the current paper, when the matrix
11
A is an invertible 3 by 3
matrix,
0 3 3
P I so that
1
P 0 , and
T
1 21 O,A
P A α 0. The orbital control force will not then
depend on the attitude of the follower satellite. More explicitly, it is obtained as
. t m
C -1
11 1 1
F A b a (3.58)
Eq. (3.58) is quite natural in that m
-1
11 1
A b is the total force to satisfy the full orbital
constraint and m
1
a is just the force resulting from the unconstrained system, so
m
-1
11 1 1
A b a is the additional control force that is necessary to satisfy the given orbital
constraints. Also, from Eq. (3.41),
1
m
-1
-1 T
22 u 22
D = A M A and the attitude control torque is
explicitly given by
1
1
1
2
1 1
,
2
t
m
-1
C T -1 T -1
22 22 u 22 2 22 2 21 11 1
-1
T -1 T C
22 22 u 22 2 22 2 21 1 21
Γ E A A M A b A a A A b
E A A M A b A a A a A F
(3.59)
77
where Eq. (3.58) has been used in the last equality and Eq. (3.59) shows the relationship
between the orbital control forces t
C
F and the attitude control torques t
C
Γ . If Eq.
(3.49) is employed, the attitude control torque is alternatively given by
2
1
2 .
t
m
-1
C T T -1 T -1
1 22 22 1 1 22 2 22 2 21 11 1
-1
T T -1 T C
1 22 22 1 1 22 2 22 2 21 1 21
Γ E A A E J E A b A a A A b
E A A E J E A b A a A a A F
(3.60)
From Eqs. (3.59) and (3.60), it is noted that the attitude control in general depends on the
follower satellite’s orbital tracking requirements.
Such a situation (
11
A is invertible) would occur, for example, if instead of Eq. (3.8a)
the two constraint equations
0
cos y t and
0
sin z t , where
0
is a given
orbital frequency in the yz-plane of the Hill frame, are used to specify the orbital
requirements for the ith follower, so that along with Eq. (3.8b) there exist three
independent orbital constraints (requirements). The invertibility of
11
A assures that the
three orbital requirements are independent of one another and the control force needed to
satisfy them then becomes independent of the follower satellite’s attitude. In brief, if
there are three independent orbital constraints, then there is no need to consider any
attitude constraints in designing the orbital controller.
Case 2: When
T
1 21
P A 0 .
First it is noted that the third column of the 3 by 3 matrix
T
1 21
P A is always zero because
T T
21 21 3 1
A A 0 . For the controller designed in the current chapter, the orbital control
78
is independent of the attitude control (or any attitude constraints) at each instant of time t
when there exists a 3 by 3 matrix H such that
.
T
21 0
A P H (3.61)
The following is the derivation of Eq. (3.61). The solution
T
21
A to the matrix equation,
T
1 21
P A 0 , is given by (Udwadia and Kalaba [44])
,
T +
21 3 3 1 1
A I P P H (3.62)
where H is an arbitrary 3 by 3 matrix. With the fact that
+
+
+ + +
1 3 3 0 3 3 11 11 3 3 11 11 1
P I P I A A I A A P and the matrix
1
P is idempotent,
the right hand side of Eq. (3.62) becomes
0
P H , which proves Eq. (3.61). With Eq. (3.61),
it follows that
T 2
1 21 1 0 3×3 0 0 0 0 0 0
P A P P H I P P H P P H P P H 0 (3.63)
so that the separation principle holds. This case is more limited than the previous one in
that Eq. (3.61) may be satisfied only at several time points during the maneuver. In
Appendix C, it is shown that the 2 by 3 matrix
21
A always has full rank, that is, rank of
2. Then, the rank of the matrix
T
21
A is 2 and the rank of the matrix
0
P is s , so it is
concluded that from Eq. (3.61)
, rank rank rank
T
21 0 0
A P H P (3.64)
which yields
2 . s (3.65)
79
Eq. (3.65) states that for a matrix H to exist the number of independent orbital
constraints should be equal to or more than two. If 3 s , the matrix
11
A at the time t is
invertible,
1
P 0 , and this falls under Case 1, and the separation principle always applies.
If 2 s , then 3 1 rank s
1
P from Eq. (3.44), so that the rank of the 3 by 3 matrix
T
1 21
P A is 0 or 1. The rank of 0 means that
T
1 21
P A is a null matrix, Eq. (3.57) is satisfied,
and the separation principle holds. The rank of 1 means that the first two nonzero
columns of the 3 by 3 matrix
T
1 21
P A are linearly dependent and the separation principle
holds only at the time instants when there exists a 3-vector
O,A
α such that
T
1 21 O,A
P A α 0.
Finally, if 1 s , a matrix H does not exist and
T
1 21
P A is never zero, and the rank of
T
1 21
P A is 1 or 2 (since the third column of
T
1 21
P A is always zero). When the rank is 1, the
separation principle holds only at the instants when
T
1 21 O,A
P A α 0 is satisfied as before.
When the rank is 2,
T
1 21 O,A
P A α 0 as long as
O,A
α 0 because no linear combination of
the two nonzero columns of
T
1 21
P A makes
T
1 21 O,A
P A α zero, which means that the
separation principle never holds. The results are summarized in Table 3-1.
Remark: One of the strengths of the analytical solution is that it is possible to estimate in
advance the magnitude or the order of the required control force and/or torque. For
example, the second component of t
C
F is
m m
T T
1 21 O,A 1 21 2 21 1 22 2 21 O
P A α P A D b A a A a A α (3.66)
80
from Eq. (3.48). With the definitions of each term in the right hand side of Eq. (3.66), the
matrix norm of D is very small if the absolute distance between the follower satellite and
the target, that is, t t t
t
ξ X X is sufficiently large. For example, provided that
the target is set as the center of the Earth, t ξ is of the order of
6
10 (in meters) if the
follower satellite is in low-Earth orbit. Then, the matrix norm of
22
A which is almost the
same as t ξ is also large, which makes the norm of D very small because D is the
inverse of the matrix multiplication of two
22
A ’s. The norms of the elements in the
parenthesis (i.e.
2
b ,
21 1
A a ,
22 2
A a , and
21 O
A α ) of the right hand side of Eq. (3.66)
are much smaller than t ξ since the norm of
2
a is of the order of the norm of
2
u ,
and the norms
1
P and
21
A are both of the order unity. It is therefore concluded that
m
T
1 21 O,A
P A α is a lot smaller than m
O
α . In many cases, hence, the first component m
O
α
dominates the second component m
T
1 21 O,A
P A α in the control force t
C
F for orbital
control.
For reference, the attitude control torques are described by
2
m
t
C T
1 22 O,A
Γ E A α from
Eq. (3.49) and it is shown in Appendix C that the rank of the matrix
T
22 1
A E has full rank
(i.e. rank of two), so t
C
Γ is never zero as long as
O,A
α is not zero. Also, although
O,A 2 21 1 22 2 21 O
α D b A a A a A α may be small as stated earlier, the torque t
C
Γ is
not negligible because the matrix norm of
22
A is large (of the order of ξ ).
81
1 s
1 rank
T
1 21
P A
The separation principle holds only at the time sections when
T
1 21 O,A
P A α 0.
2 rank
T
1 21
P A The separation principle never holds as long as
O,A
α 0.
2 s
0 rank
T
1 21
P A The separation principle holds.
1 rank
T
1 21
P A
The separation principle holds only at the time points when
T
1 21 O,A
P A α 0.
3 s
0 rank
T
1 21
P A The separation principle always holds.
Table 3-1. Various conditions for the separation principle depending on the rank s of
11
A
3.3 Numerical Example
In this section, an example is introduced to demonstrate the applicability of the control
methodology suggested in the previous sections. It is straightforward to extend this
example for applications to more general situations. The numerical integration is done in
the Matlab environment, using a variable time step (4,5)-modified Runge-Kutta
integrator with a relative error tolerance of
13
10
and an absolute error tolerance of
16
10
.
Let us consider a formation system in which there is one follower satellite whose
mass is 6.5 kg m . Also, its moments of inertia along its respective body-fixed axes are
taken to be
2
1.88 1.88 1.1 kg m diag J . The nonuniform Earth gravity
perturbation up to the 4
th
degree (i.e. 4 l in Eq. (2.53)) is assumed to apply to the leader
82
and follower satellites, and the initial radius of the leader satellite’s orbit is
6
7.1781363 10 m
L
r (i.e. the altitude is 800 km ). The initial right ascension of the
ascending node (
L
) and the inclination (
L
i ) of the leader satellite’s orbit are 30
and
80
, respectively, and it is assumed that the leader satellite is at the ascending node at the
initial time ( 0 t ). The leader satellite’s ‘nominal’ orbital frequency and its ‘nominal’
orbital period (the orbit is not closed) are given by, respectively,
3 3
3
2
1.0381292 10 rad/s, 6.0524127 10 s 1.6812257 hr,
L L
L L
GM
n P
r n
(3.67)
and 2
L
P is chosen as the duration of time used for numerical integration. Figure 3-1
shows the leader satellite’s trajectory in the ECI frame. The trajectory in 3-dimensional
space is given in the upper left figure. The upper right, the lower left, and the lower right
figures show the trajectories projected onto the XY -, YZ -, and XZ -planes, respectively.
The trajectory is not closed due to the nonuniform Earth gravity perturbation although it
is difficult to see this from the figure.
83
Fig. 3-1 The leader satellite’s trajectory in the ECI frame
For orbital control the follower satellite which is also under the nonuniform Earth
gravity perturbation up to the 4
th
degree is required to satisfy Eqs. (3.8a) and (3.8b), or
stay on a circle with a constant radius 2 km when projected onto the yz plane in
the Hill frame (centered at the leader satellite). As an attitude requirement the z axis of
the body frame of the follower satellite, is required to point to the center of the leader
satellite at all times, that is, t
t
X is the position vector of the leader satellite which is
84
time-varying. One of the correct initial conditions for orbital motion of the follower
satellite in the Hill frame is given by
0 0 m, 0 2000 m, 0 0 m, x y z
0 1.0381290 m/s, 0 0.02076223 m/s, 0 2.0762581 m/s. x y z (3.68)
However, since these relative initial position and velocity should be measured, for
example, from a carrier-phase differential GPS sensor, they may be corrupted by sensor
noise. For example, with filtered carrier-phase differential GPS signals, the relative
position measurements are predicted to be accurate to approximately 2-5 cm and the
velocity noise is predicted to be on the order of 2-3 mm/s (Busse et al. [45]). In the
current chapter, for illustration purposes let us consider measurement errors in excess of
these estimates. The errors in the initial position are taken to be 6 cm in the radial
direction (x direction in the Hill frame) and 6 cm in the intrack direction (y direction
in the Hill frame), and the errors in the initial velocity are taken to be 2 mm/s and
2 mm/s , respectively, in the radial and intrack directions, that is,
0 0.0244949 m, 0 2000.0244949 m, 0 0 m, x y z
0 1.0401290 m/s, 0 0.002 m/s, 0 2.0762581 m/s, x y z (3.69)
and Eq. (3.69) can be transformed into the ones in the ECI frame using Eq. (2.30):
6 6 3
0 6.2162748 10 m, 0 3.5893689 10 m, 0 1.9696396 10 m, X Y Z
0 646.8736516 m/s, 0 1118.3457218 m/s, 0 7338.9841623 m/s. X Y Z
(3.70)
85
The initial conditions for attitude motion are chosen as
0 1 2 3
0 0.07073720, 0 0.9915440, 0 0.1024636, 0 0.0365786, u u u u
4 6 4
0 1 2 3
0 5.6671227 10 , 0 2.8402648 10 , 0 3.6375241 10 , 0 0. u u u u
(3.71)
These initial conditions do not satisfy the attitude constraints Eqs. (3.19c) and (3.19d),
but do satisfy the unit-norm constraint Eq. (3.30) to describe a physical rotational motion.
Thus sensor errors in attitude measurements are also included. More specifically, if is
defined as the angle between the vector ˆ
body
z and the vector Tξ connecting the follower
satellite and the target (the leader satellite), that is,
ˆ : ,
body
z Tξ (see Fig. C-1 in
Appendix C), then is approximately 0 1.02
at the initial time; exact satisfaction
of the attitude constraints would require to be zero (see Eqs. (3.19)). It should be noted
that since sensor measurement noise is present affecting the correct values of the initial
conditions, both in orbit and attitude, Eq. (3.46) should be used to obtain the control force
and torque (instead of Eq. (3.52)) and Eq. (3.53) should be used to obtain the controlled
position and velocity history of the follower satellite (instead of Eq. (3.54)). The
parameter values, 0.002
i
and 0.002
i
( 1,2,3,4,5 i ) are chosen in the constraint
equations Eq. (2.48).
Figure 3-2 represents the orbits of the follower satellite with control projected on the
yz plane and xz plane in the Hill frame, respectively. The scale is normalized by
0
and these trajectories are obtained by integrating twice the analytical results given in Eq.
86
(3.53). Although it is not obvious from this figure, there are small initial discrepancies
from the desired formation configuration of the follower satellite due to the incorrect
initial conditions caused by sensor measurement errors. Despite these measurement
errors, asymptotic convergence to the desired orbital and attitude trajectories occurs as
seen in Fig. 3-8, which will be discussed later. As seen in Fig. 3-2, the follower satellite
quickly moves on the constraint circle of (normalized) radius unity and satisfies the linear
constraint as required by Eqs. (3.8). Figure 3-3 depicts the orbits of the follower satellite
without control projected on the yz plane and xz plane in the Hill frame, respectively.
As seen from the figures, without control, the follower satellite appears to show
unbounded motion and is moving leftward in the y direction, necessitating control.
In Fig. 3-4, the time histories of each component of the quaternion for the follower
satellite with control are shown where time is normalized by
L
P , the nominal period of
the leader satellite (see Eq. (3.67)). Figure 3-5 depicts quaternion time histories for the
follower satellite without control. As expected, each component starts with its initial
quaternion value and with the initial slope of its quaternion derivative, and maintains this
oscillation in the absence of any control force and/or torque.
87
Fig. 3-2 Controlled orbital motion in the yz- and xz-planes of the follower satellite. The follower
satellite stays on the circle of radius around the leader satellite (as seen in the Hill Frame).
Fig. 3-3 Uncontrolled orbital motion in the yz- and xz-planes of the follower satellite (as seen in the
Hill Frame). The initial and terminal points indicate the start and end of the integration interval.
88
Fig. 3-4 Quaternion of the follower satellite with control
Fig. 3-5 Quaternion of the follower satellite without control
Figure 3-6 depicts the control forces t
C
F per unit mass of the follower satellite in
order to maintain the desired orbital formation. These forces are directly obtained using
the analytical results, Eq. (3.46). The force components are described in the ECI frame,
89
and time is normalized by the nominal period of the leader satellite (i.e.
L
P ). As seen,
relatively large control forces are brought into play initially in order to eliminate the
initial errors in the orbit of the follower satellite. After this initial adjustment, as seen
from the figure, the forces required for orbital control become small. As noted before,
one can reduce the magnitude of the control force by choosing smaller values of
i
’s
and
i
’s than those used here; but this would, of course, come at the expense of a longer
time duration needed to get the follower satellite to its required trajectory. Figure 3-7
illustrates the control torques t
C
Γ per unit mass of the follower satellite for satisfying
the attitude requirements. These torques are also directly obtained using the new
analytical results, Eq. (3.46). The torque components are described in the body frame of
the follower satellite, and it is seen again that relatively large control torques are required
to mitigate initial attitude errors. From the figure, little torque along the z-axis is needed
once the target pointing becomes accurate. This means that only the x and y components
are necessary to control the orientation of the follower satellite. This result can be derived
from the analytical solution given in Eq. (3.49). From Eq. (C-4) of Appendix C, it is seen
that the third row of the matrix
T
1 22
E A is zero because 0
1 2
Tξ T ξ if the attitude
constraints are satisfied by Eqs. (3.19c) and (3.19d), so the z-component of the control
torque is identically zero.
90
Figure 3-8 represents errors in the satisfaction of the desired trajectory requirements
in the presence of sensor measurement errors. The first error plotted is
2 2
1
: e t
2 3
R X R X (see Eq. (3.9a)) and the second error
2
: 2
L
e t r
1 3
R X R X (see Eq. (3.9b)), both being measured in meters. The third
error is the pointing error. As previously defined, is the angle between the vector ˆ
body
z
and the vector Tξ so that
3
: e t in degrees, which is required to be zero. The fourth
error is the unit-norm quaternion constraint,
2 2 2 2
4 0 1 2 3
1 e t u u u u . Initially the
first three errors are large since the incorrect initial conditions, because of sensor
measurement noise, are used; however, they asymptotically converge to zero, as
expected, as time progresses. The small boxes that magnify the errors over the duration
1.95 2
L L
P P confirm this; the final errors in the satisfaction of the desired trajectory
requirements are of the order of
7
10
,
7
10
,
6
10
, and
15
10
, respectively. Even in the
presence of the incorrect initial conditions, the final orbital tracking errors are remarkably
small, when compared with existing results such as shown in How and Tillerson [46] in
which linear programming is employed and the follower satellite is required to stay
around the leader satellite within an errors box with dimensions of several meters, even
considering the differences between control methodologies and formation configuration
requirements. Also, for comparison, the same errors are plotted in Fig. 3-9 in which no
91
control forces and no control torques are applied. As expected, it is seen that the errors
are very large by several orders of magnitude.
The analytical results obtained in the previous section related to the separation
principle can provide new insights into the control dynamics. Since there are two
independent orbital constraints Eqs. (3.8a) and (3.8b) in this example (i.e. 2 s ), from
Table 3-1 the rank of
T
1 21
P A should be 0 or 1 during the maneuver. When the rank is zero,
the separation principle holds and when the rank is one, the two nonzero column vectors
of
T
1 21
P A are linearly dependent and the separation principle holds only at the instants
when
T
1 21 O,A
P A α 0 is satisfied. According to the numerical computation, the rank of
T
1 21
P A is always 1 throughout the maneuver in this example, and Fig. 3-10 shows that the
separation principle appears to apply only at isolated times, as happens for the first time
at 0.0058367
L
t P , at which the three components of the vector
T
1 21 O,A
P A α (the second
term on the right hand side in Eq. (3.38)) simultaneously become zero. At the beginning
the separation principle applies relatively frequently (see the upper inset box in Fig. 3-10)
but as the errors become small, the frequency with which it applies decreases (see the
lower inset box). At times other than when
T
1 21 O,A
P A α 0 the orbital control is not,
strictly speaking, independent of the attitude of the follower satellite, though the
contribution of
T
1 21 O,A
P A α is several orders of magnitude smaller than that of the first
92
term
O
α (see Eq. (3.46)). Comparison with Fig. 3-6 shows that the first term,
O
α , which
is of the order of
4
10
, mostly determines the orbital control force t
C
F .
Fig. 3-6 Required control forces per unit mass in the ECI frame for the follower satellite from
Eq. (3.46)
Fig. 3-7 Required control torques per unit mass in the body frame for the follower satellite from
Eq. (3.46)
93
Fig. 3-8 Errors in the satisfaction of the orbital and attitude requirements for the follower satellite
with control in the presence of sensor measurement errors showing asymptotic convergence to the
desired orbital and attitude requirements.
Fig. 3-9 Errors in the orbital and attitude requirements when no control forces and no control
torques are applied to the follower satellite.
94
Fig. 3-10 Each component of the
T
1 21 O,A
P A α to check if the separation principle holds
As a final application of the approach developed in this chapter the question of
robustness under unmodelled dynamics is addressed. Though a full treatment of this topic
would take us far afield, here only ‘small’ unmodelled dynamics is considered and a
simple approach is presented for robust control. In the next chapter a more advanced
controller shall be designed, assuming unknown, but bounded, model uncertainties.
The analytical results obtained herein can then be thought of as reference inputs in the
presence of small external disturbances. In the presence of disturbances, the displacement
vector t q of the follower satellite will differ from its reference displacement t
ref
q ,
which is obtained by integrating Eq. (3.53) twice assuming no disturbances. If the error,
: t t t
ref
e q q , is small, then linear control theory would apply, and the closed-
95
form analytically obtained control force (obtained from Eq. (3.46)) can be simply
augmented by a PD controller so that
,
D P
K K
C
Mq Q Q Q e e
(3.72)
where M is the mass matrix, q is the generalized displacement vector which is required
to track
ref
q , Q is the given (impressed) force,
C
Q is the nominal control force and
torque analytically obtained by Eq. (3.46), and Q
is the unmodelled external disturbance.
The parameters
D
K and
P
K are the control gains of the PD controller. Figure 3-11
illustrates a block diagram showing the error correction that is embedded in Eq. (3.72).
For example, assume that the following orbital disturbance force Q
applies to the
follower satellite:
0.1 sin 0.1 cos 0 0 0 0 ,
D D L D D L D
f f n t f f n t f
T
Q
(3.73)
where
6
2.418 10 N
D
f
and
L
n is the ‘nominal’ orbital frequency of the leader
satellite given in Eq. (3.67). Figure 3-12 shows the results when using an add-on PD
controller (Eq. (3.72)) with gains of 0.02 kg/s
D
K and
2
0.5 kg/s
P
K to compensate
for the disturbance Q
given by Eq. (3.73). As seen, each component of the position
vector is bounded with steady-state tracking errors of the order of
5
10
meter. Figure 3-
13 shows the control forces per unit mass exerted by the PD controller. The magnitude of
each component is much smaller than that shown in Fig. 3-6, and hence adding this PD
controller to the original controller designed using Eq. (3.46) would lead to robust control
96
provided, of course, that the control design requirements permit the presence of such
small steady-state tracking errors.
Fig. 3-11 Block diagram showing compensation for unmodelled dynamics using a PD controller for
robust control
Fig. 3-12 Tracking errors in position showing robust control in the presence of small unmodelled
dynamics. Note the vertical scale.
97
Fig. 3-13 Control forces per unit mass exerted by the PD controller
98
Chapter 4
METHODOLOGY FOR SATELLITE FORMATION-
KEEPING IN THE PRESENCE OF SYSTEM
UNCERTAINTIES
4.1 Introduction
In this chapter the formation-keeping control methodology in the presence of model
uncertainties is handled. It was shown in Chapter 3 that under a small external
disturbance, we could combine the obtained solutions with a simple PD controller to
mitigate the effect of this disturbance. However, in this case the magnitude of the
disturbance was small and known, and there still existed tiny steady-state errors. Instead,
it is assumed in this chapter that there are model uncertainties which are unknown, but
bounded, and it is shown that the nonlinear controller from the previous chapters is
augmented by an additional additive controller based on a generalization of the notion of
sliding surface control. This then provides a general approach to the control of the
uncertain satellite system, leading to a set of closed-form nonlinear controllers that
satisfies the desired attitude and orbital requirements within desired error bounds.
99
4.2 The Description of Constrained Mechanical Systems
As stated earlier, we denote the nominal system as our best assessment of the actual real-
life system, that is, it is assumed that no uncertainty is involved in the system. It is useful
to conceptualize the description of such a nominal multi-body system in a three-step
procedure, which is similar to the one handled in the previous chapters. We do this in the
following way:
First, we describe the so-called uncontrolled (unconstrained) system in which the
coordinates are all assumed to be independent of each other. The equation of motion of
this system is given, using Lagrange equation, by
, , , , t t M q q Q q q (4.1)
with the initial conditions
0 , 0 , t t
0 0
q q q q
(4.2)
where q is the generalized coordinate n-vector and t is time; 0 M is the n by n mass
matrix which is a function of q and t, and Q is an n-vector, called the ‘given’ force,
which is a known function of q, q , and t . From Eq. (4.1) we find that the acceleration of
the uncontrolled system is given by
: , , , . t t
-1
a M q Q q q
(4.3)
100
Second, we impose a set of control requirements as constraints on this uncontrolled
system. We suppose that the uncontrolled system is now subjected to the p constraints (or
trajectory/orientation requirements) given by
( , , ) 0, 1,2, , ,
i
t i p q q (4.4)
where ( ) r p equations among
1
,
2
, …,
p
in Eq. (4.4) are functionally independent.
The constraints described by Eq. (4.4) include all the usual varieties of holonomic and/or
non-holonomic constraints. We shall assume that the initial conditions Eq. (4.2) satisfy
these p constraints. Therefore, the components of the n-vectors
0
q and
0
q
cannot all be
independently assigned. We further assume that the set of trajectory requirements given
by Eq. (4.4) is smooth enough so that we can differentiate them with respect to time t to
obtain the relation
, , , , , t t A q q q b q q (4.5)
where A is a p by n matrix whose rank is r, and b is a p-vector. We note that each row
of A arises by appropriately differentiating one of the p constraint equations in the set
given in relation Eq. (4.4).
Using the information in the previous two steps, in the last step we obtain the
description of motion of the ‘controlled nominal system,’ or the ‘nominal system,’ for
short as
101
, , , , , , t t t
c
M q q Q q q Q q q
(4.6)
where
c
Q is the control force n-vector that arises to ensure that the control requirements
Eq. (4.5) are satisfied. The explicit equation of motion of the nominal system is given by
the fundamental equation
( ) ( ),
T -1 T +
Mq Q A AM A b Aa (4.7)
wherein the various quantities have been defined in the previous two steps and the
superscript “+” denotes the Moore-Penrose (MP) inverse of a matrix. In the above
equation, and in what follows, we shall suppress the arguments of the various quantities
unless required for clarity.
The control force that the uncontrolled system is subjected to, because of the presence
of the control requirements Eq. (4.4), can be explicitly expressed as
: ( ), ( ), ( ) ( ). t t t t
c c T -1 T +
Q Q q q A AM A b Aa
(4.8)
We refer to the system described by Eq. (4.7) as the ‘nominal system,’ implying that it
includes our best assessment of the information we have regarding the model of the
(follower) satellites, and we use the closed-form control force explicitly given by Eq.
(4.8) to control it so that it satisfies the prescribed trajectory/orientation requirements
given by Eq. (4.4). Pre-multiplying both sides of Eq. (4.7) by
-1
M , the acceleration of the
nominal system that satisfies the constraints Eq. (4.4) can be expressed as
( ) ( ) : , t
-1 T -1 T + -1 c
q a M A AM A b Aa a M Q
(4.9)
102
a relation which we shall require later on.
The generalized control force given in Eq. (4.8) is predicated on perfect knowledge of
the system and the use of an accurate model. Since in real-life situations uncertainties
exist, this control force needs to be modified to compensate for these uncertainties.
As an example, consider uncertainties in the modeling of the masses and moments of
inertia of the follower satellites. Such uncertainties in the mass could result from
modeling errors in fuel consumption especially for long-duration missions. Errors in
assessing the mass distribution may result in uncertainties in the determination of the
moments of inertia. The latter could substantially affect the attitude control needed for,
say, achieving the desired target pointing accuracy. These uncertainties affect the
elements of the n by n matrix M . Since the mass matrix
,t M q , in general, is a
function of the coordinate q and time t, the changes in the mass and moments of inertia
of the system cause changes in the coordinate q at each instance of time. Also, since the
given force
, ,t Q q q is also a function of the coordinate (and its derivative), uncertainty
in our knowledge of the mass and moments of inertia will, in general, proliferate into the
given force acting on the system.
Thus, in order to ensure that the follower satellites, whose models are not exactly
known, track the trajectory requirements of the nominal system, that is, the requirements
of our best-estimate system, Eq. (4.7) has to be replaced with
103
, , , , t t t t
c u
a c c a c c
M q q Q q q Q Q (4.10)
where
c
q is the generalized coordinate n-vector of the controlled actual system and
u
Q is
the additional control force n-vector that compensates for the fact that the model is
known only imprecisely, which we shall develop in closed form. The n by n matrix
: 0
a
M M δM is the actual mass matrix of the real-life system which is a function of
c
q and t, δM is the uncertainty in the mass matrix which may include, among others,
say, uncertainties in the masses and moments of inertia of the satellites; the actual ‘given’
force vector is taken to be :
a
Q Q δQ where the n-vector Q denotes the nominal
‘given’ forces, and δQ denotes the n-vector of the changes in the ‘given’ force that are
caused by the presence of the uncertainties; and we denote the unconstrained acceleration
of the actual uncertain system as :
-1
a a a
a M Q .
We now refer to Eq. (4.10) as the description of the ‘controlled actual system,’ or
‘controlled system,’ for short, implying that in addition to the control force ( ) t
c
Q
given
by Eq. (4.8) and obtained on the basis of the corresponding nominal system, the system is
also subjected to the additional control force
u
Q
to compensate for the uncertainties
involved. Pre-multiplying both sides of Eq. (4.10) by
-1
a
M , the acceleration of the
controlled system can be expressed as
. t
-1 c -1
c a a a c
q a M Q M Mu (4.11)
104
Here :
-1
a a a
a M Q and :
u
c
Q Mu , where
c
u
is the additional generalized acceleration
provided by the additional control forces
u
Q to compensate for uncertainties in our
knowledge of the actual system, which shall be developed later.
4.3 Formation-Keeping Equations of Motion: The Controlled Nominal
System
It is assumed that there are N follower satellites and that either a real or a fictitious
leader satellite leads this N -satellite formation. The i
th
follower satellite has a nominal
mass
( ) i
m and has a diagonal inertia matrix,
( ) i
J , in which the nominal moments of
inertia along its body-fixed principal axes of inertia are placed.
As before, it is assumed that the position vector of the center of mass of the i
th
follower satellite in the Hill frame is given by
( ) ( ) ( ) ( ) i i i i
x y z
T
x
and its orientation
is described by the quaternion
( ) ( ) ( ) ( ) ( )
0 1 2 3
i i i i i
u u u u
T
u
. Then we define the
generalized displacement 7-vector as
( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )
0 1 2 3
, 1,2, , .
i i i i i i i i
t x y z u u u u i N
T
q
(4.12)
If we assume there is no perturbation applied to both the leader and follower
satellites, then we have the following equation of uncontrolled motion for each follower
satellite from Eq. (2.73):
105
( ) ( ) ( )
( )
( ) ( ) ( ) ( ) ( ) ( ) ( )
,
1
2
i i i
i
Hill
i i i i i i i
t
N
T -1
-1
1
a M Q
a
E J ω J ω u u
(4.13)
where the 7 by 7 mass matrix is
( )
( )
( )
,
i
i
i
m
3 3 3 4
4 3 u
I 0
M
0 M
(4.14)
and
( ) ( ) ( ) ( )
ˆ
: 4
i i i i
T
u
M E J E
is previously defined.
Eq. (4.13) describes the orbital and attitude dynamics of the i
th
follower satellite when
no control forces are applied to the system yet so that it satisfies our desired attitude and
orbital requirements.
When these trajectory/orientation requirements, which are of the form of Eq. (2.6),
are imposed, the generalized control force required to follow them is explicitly given by
Eq. (2.8). In addition, we can relate the 4 by 1 generalized quaternion torque,
( ) i
u
Γ , which
is determined by Eq. (2.8), to the 3 by 1 physically applied torque
( ) ( ) ( ) ( ) i i i i
x y z
T
Γ
, about the body axis of the i
th
follower satellite, through the
relation Eq. (2.74) or Eq. (3.6).
106
4.4 Formation-Keeping Equations of Motion: The Controlled Actual
System
Real-life multi-satellite systems usually have uncertainties, which, in general, arise due to
our lack of precise knowledge of the physical system and/or of the ‘given’ forces applied
to the system. With the conceptualization of the nominal system given in the previous
chapters, these uncertainties are now assumed to be encapsulated in the elements of the n
by n mass matrix M and the n-vector Q (see Eq. (2.3)). These uncertainties cause an
error in satisfying our desired control requirements Eq. (2.5), and a corresponding
difference between the trajectories of the real-life uncertain system and the nominal
system.
We start to define the tracking error signal as
. t t t
c
e q q (4.15)
Differentiating Eq. (4.15) twice with respect to time, we get
,
c
e q q
(4.16)
which upon use of (4.9) and (4.11) yields
, , , , , ,
: ( ) : .
t t t t t
-1 -1 c -1
a c c a c a c
-1 -1
a c a c c c
e a q q a q q M q M q Q M Mu
δq M Mu δq I I M M u δq u Mu
(4.17)
In the above equation, we have defined
107
, , , , ,
, , , , ,
t t t t t
t t t t
-1
-1
a c c c
-1
-1 -1
c c
M I M q M q I M q δM q M q
I M q M q M q δM q
(4.18)
and denoted the acceleration δq as
, , , , , , , , , , , t t t t t t
-1 -1 c
c c a c c a c
δq q q q q a q q a q q M q M q Q (4.19)
where :
-1
a a a
a M Q , with : , , , t t
a c c
M M q δM q and : , , , , . t t
a c c c c
Q Q q q δQ q q
The aim in this section is to find a suitable bound on δq which we shall use in the
following section to develop a set of additive controllers to compensate for the
uncertainties involved in our knowledge of the actual satellite system.
Using Taylor’s expansion, Eq. (4.19) can be expanded as
, ,
, , , ,
1 1
1
, ,
, , ,
1
, , , , , , , , , ,
, ( ) ( )
( ) , , (
n n
a i a i
t cj j t cj j
j j cj cj
n
a ik a i
t cj j t
j
cj cj
t t t t t
Q Q
t q q q q
q q
M Q
q q t q
q q
-1 -1
c c a a
-1
a q q q q
q a q q
δq q q q q M q Q q q M q Q q q
M q
Q q q
,
, ,
1 1
1
,
,
1
) ( )
, ( ) ,
. . ., for 1 ,..., and 1 ,..., ,
n n
a i
cj j t cj j
j j
cj
n
a ik
t cj j
j cj
Q
q q q
q
M
t q q t t
q
HOT i n k n
q q
-1 -1 c
a q
M q M q Q
(4.20)
where H.O.T. denotes the higher order terms of
c
q q and
c
q q
.
We note that in Eq. (4.20),
, a i
Q ,
cj
q , and
j
q denote the corresponding i
th
and j
th
components of the n-vectors
a
Q ,
c
q , and q , respectively. Also
1
, a ik
M
represents the (i,k)
element of the n by n matrix
-1
a
M .
108
The aim is to develop a controller
c
u
such that the motion of the controlled actual
system closely tracks the motion of the nominal system, and thereby satisfies the control
requirements Eq. (4.4). We assume for the moment that the compensating control
acceleration
c
u is capable of this and causes the trajectory of the controlled actual system
,
c c
q q
to sufficiently approximate that of the nominal system so that , ,
c c
q q q q .
Under this assumption, we take the lowest order terms in Eq. (4.20) and approximate
δq
as
, , , , , , , , , , . t t t t t t t t
-1 -1 -1 -1 c
a a a
δq q q M q Q q q M q Q q q M q M q Q (4.21)
Similarly, we take the lowest order terms in Eq. (4.18) and assume that 1
-1
M δM .
Thus by Taylor’s expansion, M in Eq. (4.18) can be approximated as
, , . t t
-1
-1 -1
M I I M q δM q M δM (4.22)
Since (Henderson and Searle [47])
, , , , t t t
-1 -1
-1 -1 -1 -1 -1
a
M q M q δM q M M I δM M δM M (4.23)
expanding Eq. (4.21) and utilizing Eq. (4.23), we obtain
,
-1 -1
-1 c
δq M δM δM M Q Q M δM δQ
(4.24)
which includes the combined effect of the uncertainties δM and δQ . By taking the norm
of Eq. (4.24), one can obtain an estimate of the bound, t , on δq
as
, t
-1 -1
-1 c
δq M δM δM M Q Q M δM δQ (4.25)
109
where t
is a positive function of time. This bound depends on δM and δQ, which in
turn depends on the state of our knowledge (or ignorance) about the real-life uncertain
satellite system. Eq. (4.25) can be used to get the general form of δq
and is also
applicable to those special situations in which either δM or δQ may be judged to be so
negligibly small as to be approximated by zero.
A further simplification of the right hand side of relation Eq. (4.25) using knowledge
on only the bounds on the uncertainties in the mass matrix, δM , and the given force,
δQ , under the assumption that 1
-1
M δM , yields
( ) 1 . t
-1 -1 -1 c
δq M δM M M Q Q δM δQ (4.26)
In the simulations that follow in Section 4.6, we obtain estimates on the bound on δq
by using relation Eq. (4.25), which seems to provide a close enough approximation to the
norm of δq obtained from Eq. (4.19).
While the framework developed here is sufficient for considering any one (or all
simultaneously) of these uncertainties, later on we consider in this chapter, for illustrative
purposes, uncertainties in the masses and moments of inertia of a follower satellite.
4.5 Generalized Sliding Surface Controller
Having obtained an estimate of the bound , our aim in this section is to develop a
compensating controller that can guarantee tracking (to within desired error bounds) of
the nominal system’s trajectory in the presence of uncertainties in the actual satellite
110
system. To do this we use a generalization of the concept of a sliding surface [48-51].
The formulation permits the use of a large class of control laws that can be adapted to the
practical limitations of the specific compensating controller being used and the extent to
which we want to compensate for the uncertainties.
Noting Eq. (4.17), the tracking error signal in acceleration can be expressed as
: .
-1
a c c c
e δq M Mu δq u Mu (4.27)
where
-1
-1
M I I M δM
has been used in the last equality. We again note that M is
unknown, since δM is unknown, and it is embedded in our controller
c
u . We shall show
that the uncertain term M will be taken care of by the proof of Lyapunov stability.
We now define a sliding surface
, t k t t s e e (4.28)
where
0 k
is an arbitrary small positive number and s is an n-vector. Our aim is to
maneuver the system to the sliding surface
ε
s Ω , whereupon by Eq. (4.28), ideally
speaking, when the size of the surface
ε
Ω is zero, we obtain the relation k e e , whose
solution exp t kt
0
e e shows that the tracking error
t e exponentially reduces to
zero along this lower dimensional surface in phase space.
Differentiating Eq. (4.28) with respect to time and using Eq. (4.27), we get
. k k
c c
s e e e δq u M u (4.29)
111
Since
c
q q
can be measured, to cancel the known term k k
c
e q q
in Eq. (4.29),
we choose the controller
c
u to be of the form
, k t t
c
u e v (4.30)
so that
. k t t s v δq M e v (4.31)
We note that t δq
. Here, we have used the bound
t that is related to the
uncertainties involved in the real-life satellite system and that is obtained from Eq. (4.25).
In what follows we shall denote to mean the infinity norm.
We shall now show that the system can indeed be maneuvered to the sliding surface
ε
s Ω when
ε
Ω is defined as an appropriately small surface around s 0
whose exact
description will be shortly discussed.
We start by considering a constant t such that
0
0
> 0,
n t
t
(4.32)
where
0 0
and 0 1 k t t n t M e M (4.33)
are any arbitrary positive constants over the time duration over which the control is
applied.
We now define a control n-vector
t v so that
112
: . t t v f s (4.34)
The i
th
component,
i
f s , of the n-vector
f s is defined as
/ , 1, . . .,
i i
f s i n s g (4.35)
where
i
s is the i
th
component of the n-vector s, is defined as any (small) positive
number and the function /
i
s g is any arbitrary monotonically increasing, continuous,
odd function of
i
s
on the interval ( , ) that satisfies
0
/ , if is outside the surface ,
i i i
t k t t
f s s t
t
ε
M e
s g Ω
(4.36)
where
ε
Ω
is defined as the surface of the n-dimensional cube around the point s 0
each of whose sides has a computable length (as shown below). We note that the right-
hand side of Eq. (4.36) is always less than unity since
0
k t t M e , and hence Eq.
(4.36) will always be satisfied when 1 f s .
Result 1: The control law
k t t k t t
c
u e v e f s (4.37)
with 0 k and t v defined in Eq. (4.34) to Eq. (4.36) will cause t
ε
s Ω .
Proof: Consider the Lyapunov function
1
.
2
V
T
s s (4.38)
Differentiating Eq. (4.38) once with respect to time, we get
113
. V
T
s s
(4.39)
Substituting Eq. (4.31) in Eq. (4.39), we have
. V t t t k t t t t t t
T T T T
s v s δq s M e s M v
(4.40)
Then using Eq. (4.34) in Eq. (4.40), we obtain
, V k
T T T T
s f s s δq s M e s Mf s
(4.41)
so that
. V k
T T T T
s f s s δq s M e s M f s
(4.42)
Then using relation t δq
and noting n
T
s s , we obtain
. V n t nk n
T
s f s s s M e s M f s
(4.43)
Since f s is an odd monotonically function of s , and s and f s have the same sign,
we have
.
T
s f s s f s (4.44)
Using Eq. (4.44), thus Eq. (4.43) becomes
1
1 .
V n n t nk
n n t nk
n n t k
n
s f s M f s M e
s M f s M e
s M f s M e
(4.45)
Since
0
0
+ n t
where
0
0 1 n M , we then have
114
0
. V n t t k
s f s M e
(4.46)
Since by Eq. (4.36),
0
: 0 t t k t t f s M e outside the surface
( ) t
ε
Ω , we have
, outside the surface , V n t t
ε
s Ω
(4.47)
so that the derivative V
is negative, and we have convergence to the closed set interior to
the region enclosed by the surface t
ε
Ω .
Note that for Eq. (4.47) to be satisfied, we require relation Eq. (4.36), namely
0
/ : ,
t k t t
t
t
M e
f s g s
(4.48)
where, as noted before, 1 t . Eq. (4.48) then yields
. t
-1
s g (4.49)
In the region in which s satisfies Eq. (4.49), the Lyapunov derivative V
is negative.
This shows us that the controller Eq. (4.37) will cause t s
to decrease until it reaches the
boundary t
ε
s Ω . Further, since 1 t , and the function * g is a monotonically
increasing function, t
ε
Ω
is enclosed in an n-dimensional cube of constant size around
the point s 0 , each of whose sides has length
2 2 1 : . L t t
-1 -1
g g (4.50)
115
This gives an estimate of the n-dimensional cubical region
ε
Ω (each of whose sides is
estimated to be of constant length ) to which trajectories of the controlled actual system
will be attracted to.
Noting the fact that
t s is bounded by / 2 L
inside the surface
ε
Ω , we now have
an estimate of the error bounds given by
and , as .
2
t t t
k
e e (4.51)
Further, under the proviso 1 t t M e for 0, t where 0, is the interval over
which the control is applied, which is something that we expect, we then have
0
2 / . L t t k t
-1
g (4.52)
For ease of implementation, one could choose the function t to be a constant by
taking it to be the upper bound,
m
, so that
m
t δq for 0, t . Then Eq. (4.52)
becomes
0
2 / .
m m
L k
-1
g (4.53)
One can then, accordingly, obtain an estimate of the error bounds by replacing in the
expressions in Eq. (4.51) by the expression on the right hand side of Eq. (4.53)
and , as .
2
L
t t L t
k
e e (4.54)
Main Result: The closed-from generalized sliding surface controller for the uncertain
system,
116
0
0
,
n t
t t k
c c
a c a c a
M q Q Q Mu Q Q M e f s (4.55)
where:
(i) the control force t
c
Q is given by Eq. (4.8) and is obtained on the basis of the
nominal system;
(ii) 0 k
is an arbitrary small positive number;
(iii) f s is any arbitrary monotonically increasing odd continuous function of s on the
interval
,
as described in Eq. (4.35) with
1 f s in
ε
Ω ;
(iv) t t δq
where t is chosen based on the estimate of t δq
from Eq.
(4.25);
(v)
0
is a small positive number that satisfies
0
0 1 ; n t M (4.56)
over the time duration over which the control is done; and
(vi) under the proviso, and the expectation, that 1 M e ,
0
is chosen such that
0
, k (4.57)
will cause the actual system to track the trajectory of the nominal system within estimated
error bounds Eq. (4.53).
117
Proof: Using Eq. (4.16) in Eq. (4.27), we have
,
-1
c a c
e q q δq M Mu (4.58)
so that
.
-1
c a c
q q δq M Mu (4.59)
Consider Eq. (4.19),
.
t
t t
t
-1 -1 c
a a
-1 c -1 c
a a
-1 c
a a
δq a a M M Q
a M Q a M Q
a M Q q
(4.60)
In the last equality above, we have used Eq. (4.9).
Substituting (4.60) in (4.59), we then get
. t
-1 c -1
c a a a c
q a M Q M Mu (4.61)
Pre-multiplying both sides of Eq. (4.61) by
a
M , we obtain
. t
c
a c a c
M q Q Q Mu (4.62)
Finally, using Result 1 (Eq. (4.37)) and Eq. (4.55), the main result follows.
4.6 Numerical Results and Simulations
In this section, an example is introduced to demonstrate the applicability of the control
methodology suggested in the previous sections. It is straightforward to extend this
118
example for applications to more general situations. The numerical integration
throughout this paper is done in the Matlab environment, using a variable time step
ode15s integrator with a relative error tolerance of
8
10
and an absolute error tolerance
of
12
10
.
We consider a system in which there is only one follower satellite whose nominal
mass is 120 kg m . Also, its nominal moments of inertia along its respective body-fixed
axes are taken to be
2
10 10 7.2 kg m diag J . The arbitrary value of
0
J is chosen as
2
15 kg m . By nominal we mean our best-estimate of these parameters for the actual,
real-life, follower satellite. It is assumed that the leader satellite is in a circular orbit
around a uniform spherical Earth and the radius of its orbit is
6
7 10 m
L
r . For the sake
of simplicity, the inclination (
L
i ) of the leader satellite’s orbit is taken to be 0
(that is,
the leader satellite is just above the equator), and it is assumed that the leader satellite is
on the X-axis of the ECI frame at the initial time ( 0 t ). The leader satellite’s mean
motion and orbital period are given by, respectively,
3 3
3
2
1.0780 10 rad/s, 5.8285 10 s 1.6190 hr.
L L
L L
GM
n P
r n
(4.63)
We choose 2
L
P (two orbital periods of the leader satellite) as the duration of time used
for numerical integration.
For the follower satellite, whose mass and moments of inertia we are uncertain of,
and which are in fact unknown, we use the following attitude and orbital requirements as
119
before: (1) The follower satellite’s orbit should be on a circle with radius
4
0
m 7.0 10
when projected onto the yz-plane of the Hill frame with the leader satellite located at the
center of the circle (i.e. the projected circular orbit (PCO)); and, (2) The follower
satellite, more specifically, the x-axis of its body frame, points to the center of the Earth
at all times. Besides these two requirements, we shall impose the additional constraint
that the quaternion 4-vector of the follower satellite should have unit norm. These
trajectory requirements, which are very similar to Eqs. (3.8), (3.19), and (3.30), are
summarized as:
0 0
2 , cos , sin , x z y t z t (4.64a)
ˆ 0,
X
Y
Z
body
x T
(4.64b)
2 2 2 2
0 1 2 3
1, N u u u u u (4.64c)
where x, y, z denote the position in the Hill frame of the follower satellite, and X, Y, Z the
position in the ECI frame of the follower satellite. Also, in Eq. (4.64a), three independent
orbital constraints are used instead of two constraints as in Eqs (3.8). The rotational
frequency
of the follower satellite about the leader satellite in the Hill frame (see Eq.
(4.64a)) is set to equal
L
n .
We choose the initial conditions for orbital motion of the follower satellite as
120
4
0 0 m, 0 7.0 10 m, 0 0 m,
0 37.7347 m/s, 0 0 m/s, 0 75.4695 m/s,
x y z
x y z
(4.65)
and the initial conditions for attitude motion as
4
0 1 2 3
4 4
0 1 2 3
0 0.0707372, 0 0.997482, 0 0.00498729, 0 3.536772 10 ,
0 0.00870185, 0 6.143960 10 , 0 5.403876 10 , 0 0.
u u u u
u u u u
(4.66)
It must be noted that the initial conditions given by Eqs. (4.65) and (4.66) satisfy the
trajectory requirements given in Eqs. (4.64). The generalized control force required to
satisfy these requirements is explicitly given by Eq. (2.8). Application of this control
force yields the trajectories of the motion of the follower satellite.
Figure 4-1 represents the orbit of the follower satellite projected on the yz -plane
(left) and xz -plane (right) in the Hill frame, respectively. This is the trajectory of our
nominal system. The scale is normalized by
0
, and as seen in the figure, in these
normalized coordinates the follower satellite stays on a circle of radius unity around the
leader satellite, which is located at the origin. In Fig. 4-2, the time history of each
component of the quaternion for the follower satellite is shown where time is normalized
by
L
P , the period of the leader satellite (see Eq. (4.63)). In Fig. 4-3 we show the obtained
control forces per unit mass of the follower satellite in order to follow the desired orbital
requirements. The force components are described in the Hill frame, and time is
normalized by the period of the leader satellite (i.e.
L
P ). Figure 4-4 illustrates the control
torques per unit mass of the follower satellite for satisfying the attitude requirements
121
(with no uncertainties). The torque components are described in the body frame of the
follower satellite, and these physically applied torques are obtained using Eq. (3.6).
Figure 4-5 represents errors in satisfying the desired nominal trajectories assuming no
uncertainties, described by Eqs. (4.64). Instead of Eq. (4.64b), we use the angle
between the x-axis of the body frame and the vector
T
X Y Z T
connecting the
follower satellite and the center of the Earth. We denote these errors by (a)
1
2 , e t x z (b)
2 0
cos , e t y t (c)
3 0
sin , e t z t (d)
4
, e t and
(e)
2 2 2 2
5 0 1 2 3
1 e t u u u u .
Fig. 4-1 Constrained motion of the nominal system with no uncertainties assumed
122
Fig. 4-2 Time history of quaternions of the nominal system with no uncertainties assumed
Fig. 4-3 Required control forces to satisfy the nominal orbital constraints
123
Fig. 4-4 Required control torques to satisfy the nominal attitude constraints
Fig. 4-5 The errors in the satisfaction with the nominal constraints
As stated earlier, an exact model of the mechanical system cannot be perfectly
obtained. To see how the response of the assumed nominal system can be altered through
the effect of the uncertainty in the modeling process, we consider for reasons of
simplicity only the uncertainties in the mass m
of the follower satellite and in its
moments of inertia
x
J ,
y
J , and
z
J . We estimate that the actual values of these
124
parameters differ from our nominal (best-estimate) values by a random uncertainty of
10% of the nominal values chosen. For illustrative purposes, we pick a specific system
with 12, m 1,
x
J 1,
y
J and 0.72
z
J and perform a simulation again using
Eq. (4.9), except that we replace our best-estimate mass matrix of the uncontrolled
system (see Eq. (4.1)) with the actual mass matrix :
a
M M δM , where M is defined
in Eq. (4.14), with all other parameter values the same as previously prescribed. We note
that the elements of the 7 by 7 symmetric matrix
a
M are given in a manner similar to Eq.
(4.14). In this case, we have replaced and
i
m J in Eq. (4.14) with m m m and
i i i
J J J respectively. Using the control forces and torque obtained under the
assumption that the mass and moments of inertia are those of the nominal system
( ) t
c
Q , we obtain
: , , . t t
-1 -1 c
a a
q M Q q q M Q
(4.67)
The trajectories ( , q q
) of the system Eq. (4.67) with the actual mass and moments of
inertia are determined.
Figure 4-6 shows these orbital trajectories of the actual system projected on the yz -
plane (left) and xz -plane (right) in the Hill frame, respectively. Figure 4-7 depicts the
time history of quaternions of the actual uncertain system. Both figures are different from
those obtained from the nominal system, showing that a miss-assessment of the mass and
moments of inertia of the follower satellite can have significant consequences. The
125
resulting quaternions satisfy neither the Earth-pointing constraint nor the unit-norm
constraint.
Fig. 4-6 Motion of the actual system when 10% uncertainties in the mass and moments of inertia
of the follower satellite are involved
Fig. 4-7 Quaternions of the actual system when 10% uncertainties in the mass and moments of
inertia of the follower satellite are involved
126
We next select the structure and parameters for the controller
c
u given by Eq. (4.37).
We choose
3
/ ,
i i
f s s (4.68)
where 0 is a suitable small number, and obtain in closed-form the additional
controller needed to compensate for uncertainties in the actual system as
3
0
0
/ .
n t
t k
c
u e s (4.68)
We note that with this choice of
3
/
i i
f s s , the region outside the surface
ε
Ω is the
region outside of the n-dimensional cube around s 0 , each of whose sides has length
1/3
0
2 /
m m
L k
(see Eq. (4.53)). In this region, Eq. (4.47) assures us
that the control given by Eq. (4.68) will cause t s
to strictly decrease, until it reaches the
boundary
ε
s Ω and remains inside this n-box thereafter.
Pre-multiplying both sides of Eq. (4.55) by
-1
a
M and using the additional controller
Eq. (4.68), we obtain the closed-form equation of motion of the controlled actual system
1/3 0
0
/
n t
t k
-1 c -1
c a a a
q a M Q M M e s (4.69)
which will cause the actual system to track the trajectory of the nominal system, thereby
compensating for the uncertainty in our knowledge of the actual system.
127
With the knowledge that there is a 10% uncertainty in the mass and moments of
inertia of the follower satellite, we use the norm of Eq. (4.24) to estimate t and
m
m
t . For our simulation we choose: 7 n ,
2
10
m
t
,
0
0.1 , 0.1 k ,
0
0.5 , and
4
10
. We note that the estimate of t is not sensitive to the
magnitude of additional control forces
u
Q in our control approach.
At the scales shown, the controlled trajectories of the follower satellite projected on
the yz -plane and xz -plane in the Hill frame fall exactly on those shown in Fig. 4-1. The
errors in tracking the attitude and the orbital trajectories of the nominal system are shown
in Fig. 4-8 and Fig. 4-9. In Fig. 4-8, we see that there are minute differences between the
position coordinates of the nominal and the controlled trajectories of the follower
satellite, and the attitude differences, as represented by the quaternion vectors, are also
small as seen in Fig. 4-9. We thus see that despite the fact that we have an uncertain
system, it can be controlled to behave in exactly the same way we would like our nominal
(best-estimate) system to behave.
Figures 4-10 and 4-11 respectively show the additional control forces and torques per
unit mass of the follower satellite in order to compensate for the uncertainties in the mass
and moments of inertia. Both additional control forces and torques are seen to be small
when compared with those obtained from the nominal system (Fig. 4-3 and Fig. 4-4,
respectively).
128
In our example, Eq. (4.53) yields
1/3
4
0
2 / 2 10 ,
m m
L k
(4.70)
so that we can calculate the error norm
3
/ 2 10 t L k
e (see Eq. (4.53)). Fig. 4-8
and Fig. 4-9 show that the errors are within the calculated error norms. We also note that
use of the specified smooth cubic function eliminates chattering.
Fig. 4-8 Orbital errors between the nominal and controlled systems
Fig. 4-9 Attitude (quaternion) errors between the nominal and controlled systems
129
Fig. 4-10 Required additional control forces to compensate for uncertainties in the mass and
moments of inertia of the follower satellite
Fig. 4-11 Required additional control torques to compensate for uncertainties in the mass and
moments of inertia of the follower satellite
130
Chapter 5
NEW SOLUTIONS TO THE INVERSE PROBLEM OF
THE CALCULUS OF VARIATIONS IN A SINGLE
DEGREE-OF-FREEDOM SYSTEM
5.1 The Solution to the Inverse Problem
From this chapter, we deal with the so-called inverse problem of finding a suitable
integrand so that upon the use of the calculus of variations, one obtains the equations of
motion for systems in which the forces are nonpotential. Finding a Lagrangian of the
given equation of motion of the mechanical system is important in areas such as the
development of approximate solutions of differential equations and various numerical
techniques.
Let , , G t x t x t
be a force acting on an unconstrained single-degree-of-freedom
(SDOF) system whose mass is taken to be unity. We assume that G is of class
1
C on
3
R . Our aim is to find the function
1 1 1
1 2
: , f t t R R R so that the functional
2
1
1 1 2 2
, , , with , ,
t
t
f t x t x t dt x t x x t x
(5.1)
when rendered stationary, yields the Euler-Lagrange equation of motion
, , . x t G t x t x t
(5.2)
131
We assume that the general solution , , x t g t a b of Eq. (5.2), passes for some
values of a and b through
1 1
x t x , and
2 2
x t x . Noting that , ,
t
x t g t a b , we
further assume that one can solve for a and b
from the expressions for x t and x t so
that , , a t x t x t
and , , b t x t x t
, where and are continuous
functions from
2 1
1 2
, t t R R . Then , ,
tt
x t g t a b , and
2 1
1 2
: , G t t R R takes
the values , , , , , , , ,
tt
G t x r g t t x r t x r
, where t , x , and r are considered to be
independent variables.
Thus we need to determine the integrand function f , i.e. , , , , t x r f t x r , such
that Eq. (5.2) is identical with the Euler-Lagrange equation,
2 2 2
2
0
f f f f
x t x t
t r x r r x
(5.3)
where the partial derivatives of f are evaluated at , , t x t x t
. Substituting for x t
for Eq. (5.2) in Eq. (5.3), we get
2 2 2
2
0
f f f f
x t G
t r x r r x
. (5.4)
Since Eq. (5.2) gives the solution for every initial condition pair
1 1
, x t x t
, the
equation
2 2 2
2
, , 0
f f f f
r G t x r
t r x r r x
(5.5)
must be an identity for all
1 1
1 2
, , , t x r t t R R . To get an explicit expression for
, , f t x r , we partially differentiate Eq. (5.5) with respect to r to obtain
132
2 2 2 2
2 2 2 2
0
f f f f G
r G
t r x r r r r r
. (5.6)
Upon defining
2
2
, , , , M t x r f t x r
r
, Eq. (5.6) becomes the partial differential
equation
0
M M M G
r G M
t x r r
(5.7)
whose solution is obtained as
, , , , , , , , , , , , , M t x r t x r t x r t t x r t x r
(5.8)
where is an arbitrary function of t, , , t x r and , , t x r , and
, , exp , , , , , ,
t
G
t a b t g t a b g t a b dt
r
(5.9)
where
, ,
, ,
, , , , , , , , 0
t
t
x g t a b
r g t a b
G
t g t a b g t a b G t x r
r r
. Finally, the function f is
found by integrating , , M t x r twice with respect to r
, , , , , , f t x r M t x r drdr r t x t x
(5.10)
where , t x and , t x are arbitrary within the requirement of Eq. (5.5). A similar
derivation for the function , , f t x r
is also given in Bolza [30], Leitmann [33], and
Leitmann [52]. The notation used here and our treatment closely follows Leitmann [52].
As a simple example, let us consider the following linear oscillator equation
2 x t cx t dx t (5.11)
133
where c and d are positive constants, with the initial conditions
0
0 x x and
0
0 x x .
For simplicity, let us assume
2
0 d c , the over-damped case. Then, it is
straightforward to show that
1 2 2 1 1 2
2 1
2 1
1
, ,
t t t t t
t x r e e e x e e r
(5.12a)
1 2 2 1 1 2
1 2 2 1
2 1
1
, ,
t t t t t
t x r e e e x e e r
(5.12b)
where
2
1
d d c and
2
2
d d c . According to Eq. (5.9), we have
2
, , .
d t
t x r e (5.13)
Using Eq. (5.10), we have the following integrand
2
, , , , , , , , , .
d t
f t x r e t x r t x r drdr r t x t x
(5.14)
Since is an arbitrary function of its arguments we may use
1 2 2 1 1 2
2 1
2 1
1
( , , ) , , .
t t t t t
t x r t x r e e e x e e r
(5.15)
Substituting Eq. (5.15) into Eq. (5.14) yields
2 1 1 2
4 2 3
2 1
2 1
1 1 1
, , , ,
2 6
t t t t d t
f t x r e e e xr e e r r t x t x
(5.16)
where , t x and , t x should be chosen so that
( )
( )
, ,
x x t
r x t
f t x r
in Eq. (5.14) satisfies
the Euler-Lagrange equation. For the sake of simplicity, we can choose the following
2 1
4 3
2 1
1 2
, 0, , .
3
t t d t
c
t x t x e e e x
(5.17)
134
Finally, the integrand function f , which yields the force G in Eq. (5.11), is given
by
2 1 1 2 2 1
2 1 1 2
4 2 3 4 3
2 1 2 1
2 1 2 1
4 2 2 3
2 1
2 1
1 1 1
, ,
2 6 3
1
3 2 .
6
t t t t t t d t d t
t t t t d t
c
f t x r e e e xr e e r e e e x
e e e x r cx e e r
(5.18)
Likewise, if , , , , t x r t x r for the linear over-damped system, one other possible
integrand function is given by
2 1 1 2 2 1
2 1 1 2
4 2 3 4 3
1 2 2 1 1 2
2 1 2 1
4 2 2 3
1 2 2 1
2 1
1 1 1
, ,
2 6 3
1
3 2 .
6
t t t t t t d t d t
t t t t d t
c
f t x r e e e xr e e r e e e x
e e e x r cx e e r
.
(5.19)
If
2
, , , , t x r t x r
is chosen, we have
2 1 2 1
2 1 2 1
2 1
2 1
2 1 2 1 2
2 2
2 2 4
2 1
6
2
3
2 1
2 1
2
6 4
2 1 2
2 1
2
2 2 2
2 1
6
2
2
4
2 1
2 1
1 1
1
2 12
, ,
1
3
4
3 2
1
12
4
t t t t
d t
t t t t
t t d t
t t
d t
t t t t t
e e x r e e r
f t x r e
e e e e xr
c
e e e x
e e x r cx
e
e e r e e e e
1
3
.
t
xr
(5.20)
We can verify that the integrands , , f t x r given in Eqs. (5.18)-(5.20) are correct by
directly substituting them into the Euler-Lagrange equation Eq. (5.3). We have thus
illustrated three (out of an infinite number of) different integrand functions f that
yield the force G given in Eq. (5.11).
135
5.2 Variational Principle for Systems with Nonpotential Forces
This section deals with the derivation of some special forms of nonpotential forces that
can be obtained through the use of the variational calculus by appropriately choosing the
integrand functions. Application of the variational calculus will thus yield the proper
equation of motion for such nonpotential systems. It is noted that in order to find
integrand functions f in Eq. (5.1) of the variational problem, the general solution
, , a t x t x t
and , , b t x t x t
must be known in order to use a prescribed
function . However, the use of , , 1 t x r , while narrowing down the available
functions f , averts this difficulty (Leitmann [33]). It is obvious then that for , , f t x r
to be explicitly determined,
, , G t x r
dt
r
in Eq. (5.9) and , , M t x r drdr
in Eq.
(5.10) should be suitably integrable, and that the main difficulties in finding the function
f
analytically therefore originate from these two integrations.
With this in mind, we next consider a general form for
, , G t x r
r
given by
1 2 3
, ,
, ,
G t x r
g t g x r g r G t x r
r
(5.21)
so that
, , G t x r
dt
r
can be easily obtained, and we can avoid the difficulty of
integrating , , M t x r drdr
by separating the variables. Eq. (5.21) leads us to assume
the following form for , , G t x r :
136
2
, , , , G t x r t r x r t x r t x (5.22)
where , , , , are arbitrary (suitably smooth) functions of their arguments. We
then explore various cases in which , , G t x r has the form given in Eq. (5.22).
As stated before, though the general solution of the nonlinear differential equation Eq.
(5.2) with , , G t x r given by Eq. (5.22) is hard to find, we can avoid the difficulty by
setting , , 1 t x r in Eq. (5.8) (Leitmann [33]). If we choose , , G t x r t r , then
, , G t x r t
r
, which is of the form of Eq. (5.21) in which
1
g t t ,
2
0 g x
, and
3
0 g r . Thus, for a nonpotential force of the form , , G t x r t r
acting on a
particle of unit mass, we can find the integrand
2
, , exp ( )
2
r
f t x r t dt
which,
when used to make the integral given in Eq. (5.1) stationary, yields the Euler-Lagrange
equation , , x t G t x t x t t x t
. In a similar manner,
2
, , G t x r x r is
another possible nonpotential force since , , 2 G t x r x r
r
has the form of Eq.
(5.21) if we let
1
0 g t ,
2
2 g x x , and
3
0 g r . As another example, for
, , , G t x r t x r , 0, r we have
1
, , , ,
d r d r
G t x r t x t x r
r dr r dr
, which is of the form of Eq.
(5.21) with
3
1 d r
g r
r dr
. This allows us to solve the inverse problem of finding
the integrand , , f t x r
which, when used to render stationary the integral given in Eq.
(5.1), yields the equation of motion , x t t x t x t
. Likewise, exploring other
137
cases such as,
2
, , G t x r t r x r , , , , G t x r t r t x r , , , , G t x r t r t x ,
2
, , , G t x r x r t x r ,
2
, , , G t x r x r t x , , , , , G t x r t x r t x ,
2
, , , G t x r t r x r t x r ,
2
, , , G t x r t r x r t x ,
, , , , G t x r t r t x r t x ,
2
, , , , G t x r x r t x r t x , and
2
, , , , G t x r t r x r t x r t x ,
a set of nonpotential forces can be found, for
which a suitable function
( )
( )
, ,
x x t
r x t
f t x r
whose integral, when made stationary, yields the
correct equation of motion. The following six sets of nonpotential forces can be adduced
from a suitable function f via the variational calculus:
2
1 1 1 1
, , , G t x t x t t x t t x t x t x t
(5.23a)
2 2 2 2
, , , , 0 G t x t x t t x t x t x t
(5.23b)
2
3 3 3 3
, , , , 0
p
G t x t x t t x t x t x t t x t x t x t
(5.23c)
2
4 4 4 4
, , ln , 0 G t x t x t t x t x t x t t x t x t x t
(5.23d)
5 5 5
, , , , 0
q x t
G t x t x t t t x t e q
(5.23e)
2
6 6 6
, , , , 0.
q x t
G t x t x t x t t x t e q
(5.23f)
Here,
i
,
i
,
i
,
i
, and
i
1,2,3,4,5,6 i are arbitrary functions of their arguments,
as indicated, and and p q are real numbers. One could relax the condition 0 x t to
0 x t in Eq. (5.23c) if p is restricted to being any real number such that
p
x t is a
real, well-defined number. Eq. (5.23a) and a special form of Eq. (5.23b), namely
2
, , , G t x t x t t x t x t
, are already known and given in Leitmann [33].
138
Furthermore, the general form given in Eq. (5.23a) includes the special examples (given
in Propositions 1-3 and Proposition 5) in Musielak [53]. We can determine , , f t x r of
Eq. (5.1) using Eq. (5.10). The resulting integrands , , f t x r of Eq. (5.1) for the forces
described in Eqs. (5.23) are respectively given by,
2
1 1 1 1
1
, , , , ,
2
f t x r t x r t x t x dx
(5.24a)
2 2
2
1
, , , f t x r drdr t x dx
r
(5.24b)
2
3 3 3 3 3 3
1
, , , , 1, 2
1 2
p
f t x r t x r t t x x dx p p
p p
(5.24c)
4 4 4 4 4 4
, , ln f t x r t x r r r t t x dx
(5.24d)
5 5 5 5 2
1
, , ,
qx
f t x r t e t t x dx
q
(5.24e)
2
6 6 6 6
1
, , erf , , 0
2
q r
e
f t x r x r q r t x x dx q
q q
(5.24f)
2
6 6 6 6
1
, , erfi , , 0
2
qr
e
f t x r x r q r t x x dx q
q q
(5.24g)
where,
1 1 1
, exp 2 t x t dt x dx
(5.25a)
3 3 3 3
exp 1 ( ) , exp 2 t p t dt x p x dx
(5.25b)
4 4 4 4
exp , exp t t dt x x dx
(5.25c)
5 5
exp t q t dt
(5.25d)
6 6
exp 2 x q x dx
(5.25e)
139
and erf and erfi denote the error function and the imaginary error function,
respectively. Eqs. (5.24f) and (5.24g) pertain to
6
G given in Eq. (5.23f). Eq. (5.24c)
does not hold when 1 p or 2 p . When 1 p , Eq. (5.23c) becomes
2
3 3 3
, , , G t x t x t t x t x t x t x t
(5.26)
and
3 3 3 3
, , ln , f t x r x r r r t x x dx
(5.27)
where
3 3
exp x x dx
. In addition, when 2 p , Eq. (5.23c) becomes
2
3 3 3
, , , G t x t x t t x t t x t x t
(5.28)
and
3 3 3 3
, , ln , f t x r t r t t x dx
(5.29)
where
3 3
exp t t dt
.
By using Eq. (5.24b) we provide in Table 5-1 some nonpotential forces
, , G t x t x t
that can be obtained ( constant ) by using the calculus of variations
along with the corresponding integrands , , f t x r . It is noted that if we substitute each
integrand
( )
( )
, ,
x x t
r x t
f t x r
into the Euler-Lagrange equation (Eq. (5.3)), we get the
corresponding force , , G t x t x t
.
140
, , G t x t x t
, , f t x r
2 2 2
x t x t
2 2 2 2 3
1
ln
3
r r r r x
,
x t
h t x t e
2
1
,
r
e h t x dx
2 2
, h t x t
x t
4 2 2
1 1
,
12 2
r r h t x dx
2 2
,
x t
h t x t e
2 2
erfi ,
2
where, erfi( ) : erf ( )
r
e
r r h t x dx
z i i z
Table 5-1. The function , , G t x t x t
of Eq. (5.2) and the corresponding function , , f t x r of
Eq. (5.10)
Other special cases such as
3
, , G t x t x t x t x t
, 0 x t ( , < 0 ) can
be obtained from the integrand
2
2
1
, ,
2
t
t
e
f t x r e x
r
by using Eq. (5.24c).
Another
integrand, using Eq. (5.24b), for the same force G is given by
1
2
tan
ln
ln 1
, , , 0
2
r
r r
r r
f t x r x r
.
As an important example, let us consider a model of the Coulomb friction force
approximated by
, , tanh , 0 G t x t x t x t x t
(5.30)
141
where 0 and 0 are constants. Comparing with Eq. (5.23b), it is obvious that
2
, t x t
and
2
tanh x t x t
. Then, Eq. (5.24b) yields the following
integrand
2 2 2 2
2 2
, , coth
1
2 ln 1 2 ln sinh , 0
2
r r
f t x r r drdr dx
Li e r r e r r x r
(5.31)
where
2
Li is the polylogarithm function, also known as the Jonquière’s function, defined
as
1
k
n n
k
z
Li z
k
(5.32)
where n is an integer. To verify that the , , f t x r in Eq. (5.31) does yield the Coulomb
friction force given by Eq. (5.30), we directly substitute it into the Euler-Lagrange
equation, Eq. (5.3). First, we have
2
2 2 2 2 2
1 2 2
1
2 2 2 ln 1 4 2 ln sinh 2 coth
2 1
r
r r
r
f re
Li e r e r r r
r e
(5.33)
where we have used the fact
1
1
n n
d
Li z Li z
dz z
. The total differentiation of Eq. (5.33)
with respect to time gives us, after some calculations
, , coth .
d
f t x t x t x t x t
dt r
(5.34)
Next, it is obvious from Eq. (5.31) that
.
f
x
(5.35)
142
With Eqs. (5.34)-(5.35), the Euler-Lagrange equation then yields the force , , G t x t x t
that is given by Eq. (5.30), when 0 x t .
143
Chapter 6
LAGRANGIANS FOR DAMPED LINEAR MULTI-
DEGREE-OF-FREEDOM SYSTEMS
6.1 Lagrangians for Damped Linear Two-Degree-of-Freedom Systems
In the present chapter, the findings obtained in Chapter 5 are extended to dissipative,
constant coefficient, linear multi-degree-of-freedom (MDOF) systems. The emphasis is
on obtaining the Lagrangians for these MDOF systems in a simple manner, using insights
obtained from our understanding of the inverse problem for SDOF and 2-DOF systems. It
is shown that the solution to a gyroscopically damped linear system is easily found as a
special case of the linearly damped case. Lagrangians for special linearly damped MDOF
systems with symmetric stiffness and damping matrices are also obtained along with the
corresponding Jacobi integrals, which are conserved over time.
We begin with the problem of finding a Lagrangian function for a linear mass-spring-
damper system in single degree of freedom whose governing equation of motion is given
by
2 0, 0, , 0 mx bx kx m b k (6.1)
where x t is a generalized displacement of the mass, the dot denotes the differentiation
with respect to time t, and m, b, and k are the mass, damping, and stiffness coefficients,
144
respectively, which are all assumed to be constants. Unfortunately, Eq. (6.1) cannot be
directly derived as the Euler-Lagrange equation from a variational principle because it
does not satisfy the Helmholtz conditions (Helmholtz [36], Santilli [37]) (see Eqs. (6.11)-
(6.14)). However, in Chapter 5 it is shown that the following Lagrangian function results
in Eq. (6.1):
2
2 2
1 1
.
2 2
b
t
m
L e mx kx
(6.2)
More precisely, substituting Eq. (6.2) into the Euler-Lagrange equation of the standard
form 0
d L L
dt x x
yields
2
2 0.
b
t
m
e mx t bx t kx t
(6.3)
Since the exponential factor in Eq. (6.3) is always positive in time, we can say that Eq.
(6.3) is ‘equivalent’ to Eq. (6.1). This exponential factor, however, plays a significant
role, because Eq. (6.3) does satisfy the Helmholtz conditions.
Thinking of the Lagrangian given by Eq. (6.2) as representing a physical, under-
damped oscillator consisting of a ‘mass’
2b
t
m
me and ‘potential energy’
2
2
1
2
b
t
m
e kx , we find
that the energy of this oscillator does not decay with time as does that of the under-
damped system described by Eq. (6.1) (Ray [54]). In fact, its energy, averaged over one
period of the oscillator, remains a constant, and hence the Lagrangian given by Eq. (6.2)
is not a physical Lagrangian, but just a convenient mathematical tool to use the
145
machinery of the calculus of variations to obtain Eq. (6.3). In a similar manner, were the
Lagangian in Eq. (6.2) to represent an over-damped oscillator with ‘mass’
2b
t
m
me
, and
‘potential energy’
2
2
1
2
b
t
m
e kx , its energy would, in general, exponentially increase with
time, contrary to behavior of the over-damped oscillator described by Eq. (6.1). It follows
then that this Lagrangian given in Eq. (6.2) does not describe the physical linearly
damped system though the two systems are described by the same equation of motion,
namely Eq. (6.1). However, a mathematical (not physical) Lagrangian may yet be useful
in areas such as the development of approximate solutions of differential equations and
various numerical techniques.
Taking a hint from Eq. (6.2), we next consider a two-degree-of-freedom (2-DOF)
mass-spring-damper system using the Lagrangian given by
2 2 2 2
1 1 1 1 2 2 2 2 1 1 2 2 1 2 1 2
1 1
,
2 2
t
L e m x k x m x k x b x x b x x dx x
(6.4)
where
i
m ,
i
k ,
i
b , 1,2 i , and and d are constants. Then, the corresponding Euler-
Lagrange equations of motion are given by
1 1 1 1 1 2 2 1 1 1 2
0, m x m x b b x k x b d x (6.5a)
2 2 2 2 2 1 1 2 2 2 1
0. m x m x b b x k x b d x (6.5b)
146
Depending on the choice of the constants, Eqs. (6.5) can represent various systems.
For example, if we choose
1
m m ,
2
2 m m , 2 ,
1 2
0 b b ,
1 2
2 k k k , and d k ,
Eq. (6.4) becomes
2 2 2 2 2
1 1 2 2 1 2
1
2 ,
2
t
L e mx kx mx kx kx x
(6.6)
and the corresponding equations of motion become
1 1 1 2
2 2 0, mx mx kx kx
(6.7a)
2 2 2 1
2 4 2 0, mx mx kx kx
(6.7b)
which describe a classically damped 2-DOF system. In fact, Eqs. (6.7) describe the
mechanical system shown in Fig. 6-1.
Fig. 6-1: Linear 2-DOF mass-spring-damper system with b = 2m
In order to get a more systematic approach to the inverse problem for a constant
coefficient linear 2-DOF system we consider the following equations of motion:
1
1 1 1 1 2 1 1 1 2
0, f x a x b x c x d x (6.8a)
147
2
2 2 1 2 2 2 1 2 2
0, f x a x b x c x d x (6.8b)
where all coefficients are constant and we have divided each equation by the
corresponding masses
1
m and
2
m . More generally, we consider the following set of
equations:
1
11 11 1 1 1 1 2 1 1 1 2
0, f x a x b x c x d x (6.9a)
2
22 22 2 2 1 2 2 2 1 2 2
0, f x a x b x c x d x (6.9b)
where 1,2
ii
i are non-zero functions of t ,
1
x ,
2
x ,
1
x , and
2
x , namely
11 11 1 2 1 2
, , , , t x x x x and
22 22 1 2 1 2
, , , , t x x x x . Next, we define functions
1,2
i
i by
1 11 1 1 1 2 1 1 1 2
: , a x b x c x d x (6.10a)
2 22 2 1 2 2 2 1 2 2
: . a x b x c x d x (6.10b)
Now let us consider the n differential equations , , , , 0
j
ij i
t q t q q q q
( , 1,2, , ) i j n , where
1 2 n
q q q
T
q is a generalized displacement n-vector
that describes the motion of a mechanical system in an n-dimensional configuration
space, and the summation convention is used for repeated indices. The question of
whether such a system can be obtained from a suitable Lagrangian, , , L t q q , through
the use of the Euler-Lagrange equations 0
i i
d L L
dt q q
1
appears to have been first
1
Since we cannot differentiate the function , , L t q q , with respect to a dependent variable, say
i
q , by
148
investigated by Helmholtz (Helmholtz [36]). The necessary and sufficient conditions for
the so-called ordered direct analytic representations are (Santilli [37])
,
ij ji
(6.11)
,
ij
ik
k j
q q
(6.12)
2 ,
j k i
ij j i k
q
q q t q
(6.13)
1
,
2
j j k i i
j i k j i
q
q q t q q q
(6.14)
where the summation convention is applied for repeated indices. Eqs. (6.11)-(6.14) are a
set of partial differential equations
2
that need to be satisfied by the 2 1 n independent
variables t, q, and q everywhere in
2 1 n
R
. Throughout this chapter we shall be dealing
with Lagrangians that provide the so-called ordered direct analytic representations of the
equations of motion (Santilli [37]).
, ,
i
L t
q
q q
we mean
, ,
i
L t
r
s q
r q
s r
, where t , s , and r are considered independent variables; similarly, by
, ,
i
L t
q
q q
, we mean
, ,
i
L t
s
s q
r q
s r
.
2
There is a slight abuse of notation here, since in Eqs. (6.11)-(6.14) the variables t, q , and q are considered to be
independent, while in the equations of motion, x and x (see Eq. (6.8)) are considered to be functions of time t.
149
Applying these conditions to Eqs. (6.9), , , 1,2 i j k , we therefore have
1
1
q x ,
2
2
q x , and
12 21
0 . Thus, Eq. (6.11) is automatically satisfied, and Eq. (6.12)
yields
11 12
2 1
0,
x x
(6.15)
21 22
2 1
0,
x x
(6.16)
from which we conclude that
11 11 1 2 1
, , , , t x x x (6.17)
22 22 1 2 2
, , , . t x x x (6.18)
Next, Eq. (6.13) yields the following three equations:
11 11 11 11
1 1 1 2 1 1 1 2 11 1 1 2
1 1 2
, 1, 1 a x b x c x d x a x x i j
x t x x
(6.19)
11 1 22 2
0, 1, 2 or 2, 1 b a i j i j (6.20)
22 22 22 22
2 1 2 2 2 1 2 2 22 2 1 2
2 1 2
. 2, 2 a x b x c x d x b x x i j
x t x x
(6.21)
Finally, Eq. (6.14) provides only one nontrivial equation for the case 1, 2 i j , or
2, 1 i j , which is
11 22
1 1 1 2 1 1 1 2 11 1 2 1 2 2 2 1 2 2 22 2
2 1
11 11 22 22
1 1 1 2 2 2
1 2
1 1
2 2
1
.
2
a x b x c x d x d a x b x c x d x c
x x
b x b a x a
t x t x
(6.22)
150
Thus Eqs. (6.17)-(6.22) must be satisfied everywhere in
5
R if Eqs. (6.9) can be obtained
from a Lagrangian ( , , ) L t x x where
1 2
x x
T
x , i.e., if
.
j
ij
i i
d L L
f
dt x x
(6.23)
We now consider two cases.
Case I: When
2
0 a
If
2
0 a , from Eq. (6.20) we have
1
22 11 1
2
, 0.
b
b
a
(6.24)
Then, from Eqs. (6.17) and (6.18),
11 11 1 2 22 22 1 2
, , , , , , t x x t x x (6.25)
and Eqs. (6.19) and (6.21) become
11 11 11
11 1 1 2
1 2
, a x x
t x x
(6.26)
22 22 22
22 2 1 2
1 2
. b x x
t x x
(6.27)
However, since
11 11 1 2
, , t x x and
22 22 1 2
, , t x x from Eq. (6.25), taking the
partial derivative of Eq. (6.26) with respect to the variable
2
x we get
11
2
0
x
. In a
similar fashion, taking partial derivatives with respect to
1
x yields
11
1
0
x
. By doing the
same with Eq. (6.27), we obtain we have
22 22
1 2
0
x x
so that
151
11 11 22 22
, , t t (6.28)
and Eqs. (6.26) and (6.27) become
11
11 1
, a
t
(6.29)
and
22
22 2
. b
t
(6.30)
The general solutions to Eqs. (6.29) and (6.30) are:
1
11 1
,
a t
t e (6.31)
2
22 2
,
b t
t e (6.32)
where
1
and
2
are non-zero constants. Substituting Eqs. (6.31) and (6.32) into Eq.
(6.24) yields
2 1 1
2 1
2
.
b t a t
b
e e
a
(6.33)
Eq. (6.33) should hold for all 0 t , so
1 2
, a b (6.34)
1
2 1
2
,
b
a
(6.35)
and
11
and
22
become
1
11 1
,
a t
t e (6.36)
1 1
22 1
2
.
a t
b
t e
a
(6.37)
Finally, we use Eq. (6.22) to get
152
2 1 1 2 1 2 1
. a d bc a a b (6.38)
In conclusion, we have the following possible case for the existence of a Lagrangian:
Case I:
1 1 1
11 1 22 1 1 2 2 1 1 2 1 2 1 2 1 1
2
, , , , 0, 0, 0.
a t a t
b
e e a b a d bc a a b a b
a
(6.39)
It is noted that Eq. (6.39) is symmetric, that is, if we interchange the coefficients of Eqs.
(6.9a) and (6.9b) with one another, Eq. (6.39) still holds.
Cases II & III: When
2
0 a
If
2
0 a , from Eq. (6.20) we have
1
0. b (6.40)
Then, Eqs. (6.19) and (6.21) become
11 11 11 11
1 2 1 1 1 1 1 2 1 11
1 2 1
, x x a x c x d x a
t x x x
(6.41)
22 22 22 22
1 2 2 2 2 1 2 2 2 22
1 2 2
, x x b x c x d x b
t x x x
(6.42)
and Eq. (6.22) is now
11 22
1 1 1 1 1 2 11 1 2 2 2 1 2 2 22 2
2 1
. a x c x d x d b x c x d x c
x x
(6.43)
However,
11 11 1 2 1
, , , t x x x and
22 22 1 2 2
, , , t x x x by Eqs. (6.17) and (6.18).
Hence, taking the partial derivative with respect to
2
x on both sides of Eq. (6.41) we find
that
11
2
0
x
. Using a similar argument for Eq. (6.42) we find that
22
1
0
x
, so that
11 11 1 1 22 22 2 2
, , , , , , t x x t x x (6.44)
153
and Eqs. (6.41) and (6.42) become
11 11 11
1 1 1 1 1 1 2 1 11
1 1
, x a x c x d x a
t x x
(6.45)
22 22 22
2 2 2 2 1 2 2 2 22
2 2
. x b x c x d x b
t x x
(6.46)
Having the conditions Eq. (6.44), we can simplify Eq. (6.43) further as
11 1 22 2
, d c (6.47)
where
11 11 1 1
, , t x x and
22 22 2 2
, , t x x from Eq. (6.44). Now let us consider two
cases.
First, if
2
0 c and
1
0 d , then the left hand side of Eq. (6.47) is a function of t ,
1
x ,
and
1
x , whereas the right hand side is a function of t ,
2
x , and
2
x . Hence,
11
and
22
can be functions only of time t , i.e.
11 11
t and
22 22
t , and Eqs. (6.45) and
(6.46) accordingly become
11
1 11
, a
t
(6.48)
22
2 22
, b
t
(6.49)
whose general solutions are respectively
1
11 1
,
a t
t e (6.50)
2
22 2
,
b t
t e (6.51)
where
1
and
2
are arbitrary non-zero integration constants. With Eqs. (6.50) and
(6.51), Eq. (6.47) now reads
154
1 2
1 1 2 2
,
a t b t
d e c e (6.52)
which is required to hold for all 0 t , so that
1
1 2 2 1 2 1
2
, , 0, 0.
d
a b c d
c
(6.53)
We thus can obtain a second possible case for the existence of a Lagrangian:
Case II:
1 1 1
11 1 22 1 2 1 1 2 2 1 1
2
, , 0, 0, , 0, 0, 0.
a t a t
d
e e a b a b c d
c
(6.54)
Finally, if
2 1
0 c d in Eq. (6.47), then from Eqs. (6.8), the equations of motion are
now
1
1 1 1 1 1
0, f x a x c x (6.55a)
2
2 2 2 2 2
0, f x b x d x (6.55b)
and the system is not coupled any more. Each of these uncoupled systems has the form
given in Eq. (6.1) for which a Lagrangian and the corresponding
ii
’s ( 1,2 i ) are well-
known and given in Eq. (6.2). Thus, this case is summarized as
Case III:
1 2
11 1 22 2 2 1 2 1 1 2
, , 0, 0, 0, 0, 0, 0,
a t b t
e e a b c d (6.56)
and the 2-DOF system is uncoupled.
Table 6-1 summarizes the three cases Eqs. (6.39), (6.54), and (6.56). Also,
corresponding to each case it includes one Lagrangian, which shall be obtained later. It is
noted that there are many other possible Lagrangians which are different from the ones
given in Table 6-1.
155
Case Conditions Lagrangian
I
1
11 1
a t
e ,
1 1
22 1
2
a t
b
e
a
,
1 2
a b ,
2 1 1 2 1 2 1
a d b c a a b ,
2
0 a ,
1
0 b ,
1
0
1
2 2 2 1 1
1 2 1
2
2 1 2 1 2
1 2 2 1 1 2
2 2
1
2 2 2
2
a t
b c
x x x
a
L e
b c b d
x x x b x x
a a
II
1
11 1
a t
e ,
1 1
22 1
2
a t
d
e
c
,
2
0 a ,
1
0 b ,
1 2
a b ,
2
0 c ,
1
0 d ,
1
0
1
2 2 2 2 1 1 1 2
1 2 1 1 1 2 2
2 2
1
2 2 2 2
a t
d c d d
L e x x x d x x x
c c
III
1
11 1
a t
e ,
2
22 2
b t
e ,
2
0 a ,
1
0 b ,
2
0 c ,
1
0 d ,
1
0 ,
2
0 (uncoupled)
1 2
2 2 2 2
1 1 1 2 2 2
1 1 1 1
2 2 2 2
a t b t
L e x c x e x d x
Table 6-1 Three cases when a Lagrangian function of the system described by Eqs. (6.8) exists and a
corresponding Lagrangian for each case
Fig. 6-2: General 2-DOF mass-spring-damper system
As a simple application let us consider a general 2-DOF mass-spring-damper system
shown in Fig. 6-2. Its equations of motion are easily obtained as
1 2 2 1 2 2
1 1 2 1 2
1 1 1 1
0,
b b b k k k
x x x x x
m m m m
(6.57)
156
2 2 3 2 2 3
2 1 2 1 2
2 2 2 2
0,
b b b k k k
x x x x x
m m m m
(6.58)
and each equation has been normalized by each mass for the correspondence with the
form given in Eqs. (6.8). The aim is to find the conditions under which a Lagrangian
would exist for this system. First, applying Case I given in Eq. (6.39) or Table 6-1, we
obtain the following conditions
1 2 2 3 1 2
2 2 2
1 2 1
, , 0.
2
b b b b b b
k b b
m m m
(6.59)
Moreover, Eq. (6.59) can be further divided into two cases whether the two masses
1
m
and
2
m are equal or not:
3 1 2 1 1 2
2 2 2 2 1 2 1 3 1 3
2 1 2 1 1
, , 0, , and , , , are arbitrary.
2
b m m b b b
b k b b m m b b k k
m m m m m
(6.60a)
3 2 1 2
1 3 2 2 2 2 1 2 2 1 3
, , 0, , and , , are arbitrary.
2 2
b b b b
b b k b b b m m m b k k
m m
(6.60b)
Next, we apply Case II given by Eq. (6.54) or Table 6-1 to obtain the conditions
3 1
2 2 1 2 1 3
1 2
0, , 0, and , , , are arbitrary.
b b
b k m m k k
m m
(6.61)
It is straightforward to check that the mechanical system shown in Fig. 6-1 falls into this
case.
Finally, we apply Case III given by Eq. (6.56) (see Table 6-1) to have the conditions
2 2 1 2 1 3 1 3
0, 0, and , , , , , are arbitrary. b k m m b b k k
(6.62)
In this case, the system is uncoupled.
157
In brief, if a given 2-DOF mass-spring-damper system described by Fig. 6-2 does not
satisfy at least one condition given by Eqs. (6.60)-(6.62), then it appears not possible to
find a Lagrangian so that the Euler-Lagrange equation 0
i i
d L L
dt x x
( 1,2 i ) yields
an ordered direct analytic representation of the equation of motion of the system. In
particular, from Eq. (6.60) we see that when the two masses
1
m and
2
m (see Fig. 6-2) are
connected solely by a damper, i.e., when
2
0 b
and
2
0 k , there appears to be no
Lagrangian
1 2 1 2
, , , , L t x x x x that provides an ordered direct analytic representation of the
equations of motion of the system described by Eqs. (6.57) and (6.58).
Up to now we have derived the necessary and sufficient conditions for which there
exists a Lagrangian function of the 2-DOF system given by Eqs. (6.8), and obtained three
cases Eqs. (6.39), (6.54), and (6.56) summarized in Table 6-1. We next address the
question of finding a Lagrangian for each of these cases. Since a Lagrangian for Case III
is already known in Chapter 5 (also in Eq. (6.2)), we focus on obtaining a Lagrangian
function for Cases I and II.
For Case I, first, from Eq. (6.39), we know that the equations of motion are given by
1 1 2
1 1 1 1 2 1 1 1 1 2
2
0,
a t
b c
e x a x b x c x a b x
a
(6.63a)
1 1
2 2 1 1 2 2 1 2 2
2
0,
a t
b
e x a x a x c x d x
a
(6.63b)
158
where
1
1 is used. We next search for the conditions under which a Lagrangian
, , L t x x exists such that the corresponding Euler-Lagrange equations yield Eqs. (6.63),
that is,
, , , , .
ij j i
i i
d L L
t x t
dt x x
x x x x
(6.64)
Expanding the total time derivative, we have the following identities (Santilli [37]):
2
,
ij
i j
L
x x
(6.65a)
2 2
,
j i
i j i i
L L L
x
x x x t x
(6.65b)
where , 1,2 i j . Knowing
ij
and
i
from Eqs. (6.39), (6.63), and (6.64), we obtain the
following Lagrangian using Eqs. (6.65):
1
2 2 1
1 2 1 1 2 2 1 2 1 2
2
1
, , , , , , ,
2 2
a t
b
L e x x x p t x x x q t x x r t x x
a
(6.66)
where
1 2
, , p t x x ,
1 2
, , q t x x , and
1 2
, , r t x x are arbitrary functions of their arguments
within the requirements that
1 1 1 1 2 1
1 1 1 1 1 2 2 1 2 2
2 1 1 2 2 2
, , .
a t a t a t
p q p r b c r q b
b e e c x a b x e c x d x
x x t x a x t a
(6.67)
For example, if we choose
1 1
2 2 1 1 2 1 2
1 2 1 2 1 2 1 2 1 1 2 2
2 2
, , , , , 0, , , ,
2 2
a t a t
c b c b d
p t x x b x e q t x x r t x x e x x x x
a a
(6.68)
then the Lagrangian in Eq. (6.66) becomes
159
1
2 2 2 2 1 1 1 2 1 2
1 2 1 1 2 2 1 1 2
2 2 2
1
,
2 2 2 2
a t
b c b c b d
L e x x x x x x b x x
a a a
(6.69)
which is shown in Table 6-1.
For Case II, following the same procedure shown in the previous Case I, we can again
obtain Lagrangian functions and one possible Lagrangian is
1
2 2 1
1 2 1 1 2 2 1 2 1 2
2
1
, , , , , , ,
2 2
a t
d
L e x x x u t x x x v t x x w t x x
c
(6.70)
where
1 2
, , u t x x ,
1 2
, , v t x x , and
1 2
, , w t x x are arbitrary functions of their arguments
within the requirements that
1 1 1 2
1 1 1 2 1 1 2
2 1 1 2 2
, , .
a t a t
u v u w v w d d
e c x d x e d x x
x x t x t x c
(6.71)
For example, if we choose
1
2 2 1 1 2
1 2 1 2 1 2 1 1 1 2 2
2
, , 0, , , 0, , , ,
2 2
a t
c d d
u t x x v t x x w t x x e x d x x x
c
(6.72)
then the Lagrangian Eq. (6.70) becomes
1
2 2 2 2 1 1 1 2
1 2 1 1 1 2 2
2 2
1
,
2 2 2 2
a t
d c d d
L e x x x d x x x
c c
(6.73)
which is listed in Table 6-1. Comparing Eqs. (6.7), which are the equations of motion of
the system shown in Fig. 6-1, with Case II, we see that
2 1 1 2 2 1 1 1 2
2
0, 0, 2, , , 1, , ,
2
k k k k
a b a b c d c d
m m m m
(6.74)
160
and the Lagrangian given in Eq. (6.73) then reduces to the one given in Eq. (6.6). The
obtained Lagrangians are summarized in Table 6-1.
6.2 Lagrangians for Special Constant Coefficient Multi-Degree-of-
Freedom Linear Systems
If the number of degrees of freedom of a system is much greater than two, then solving
the Helmholtz conditions becomes quite complex, and hence obtaining Lagrangians for
ordered direct analytic representations of n-degree-of-freedom (n-DOF) systems by
solving these conditions becomes, in general, extremely difficult, if not nearly
impossible. In this section we therefore use ideas from SDOF and 2-DOF systems to
expand our thinking to MDOF systems and hence all together bypass the need for solving
the Helmholtz conditions. For two-degree-of-freedom systems, we used the Lagrangian
given by Eq. (6.4), which we now generalize for an MDOF system as
1 1
,
2 2
t
L e
T T T
x Mx x Kx x Bx (6.75)
where
1 2 n
x x x
T
x is a generalized displacement n-vector, and M and K are
the n by n symmetric mass and stiffness matrices, respectively. Also, the scalar and the
n by n matrix B are determined according to the requirements of the problem as shall be
shown shortly. Then the equation of motion obtained via the Euler-Lagrange equation is
.
T T
Mx M B B x K B x 0 (6.76)
161
First, if we choose 0 and
T
B B (skew-symmetric), Eq. (6.76) becomes
2 , Mx Bx Kx 0 (6.77)
which is the general equation of motion for an n-DOF linear mechanical system with
gyroscopic damping. The corresponding Lagrangian is
1 1 1 1
,
2 2 2 2
L
T T T T T T
x Mx x Kx x Bx x Mx x Kx x Bx (6.78)
where M and K are symmetric matrices, and B is a skew-symmetric matrix. We note that
since 0 , the Lagrangian is a physical Lagrangian (see remarks that follow Eq. (6.3)).
As stated before, Eq. (6.77) has a skew-symmetric damping matrix, B. However, in many
practical applications, and especially in the theory of linear vibrations, the equations of
motion have symmetric stiffness and damping matrices. In order to encompass such
systems, we choose 0 and
T
B B . Then, with the Lagrangian given in Eq. (6.75),
the Euler-Lagrange equation yields
2 . Mx M B x B K x 0 (6.79)
When the skew-symmetric matrix B 0 , Eq. (6.79) reduces to the equation
, Mx Mx Kx 0 (6.80)
which describes a proportionally damped system. The Lagrangian from which this
equation is obtainable is simply given, using Eq. (6.75), by
1 1
,
2 2
t
L e
T T
x Mx x Kx (6.81)
162
for any matrix M 0 , and any symmetric matrix K . Having disposed off the case B 0
, from here on we shall concentrate then on the case when the skew-symmetric matrix
B 0 .
We would thus want the matrices 2 M B and B K in Eq. (6.79) to be
symmetric, where B 0 is a skew-symmetric matrix. The required conditions are not
obvious, so let us consider a 3-DOF system, which, by extension, will help us to adduce
the general procedure for handling linearly damped MDOF systems. We start by
considering diagonal mass matrices. If we have
1 1 12 13
2 2 12 23
3 3 13 23
0 0 0 0 0
0 0 , 0 0 , 0 ,
0 0 0 0 0
m k b b
m k b b
m k b b
M K B (6.82)
then, Eq. (6.79) becomes
1 1 1 12 13 1 1 12 13 1
2 2 12 2 23 2 12 2 23 2
3 3 13 23 3 3 13 23 3 3
0 0 2 2 0
0 0 2 2 0 ,
0 0 2 2 0
m x m b b x k b b x
m x b m b x b k b x
m x b b m x b b k x
(6.83)
and the damping and stiffness matrices of this system are not symmetric. However,
noting the negative sign in the term that involves
2
2
x in Eq. (6.69), if we choose to use the
same B matrix as in Eq. (6.82) and
1 1
2 2
3 3
0 0 0 0
0 0 , 0 0 ,
0 0 0 0
m k
m k
m k
M K (6.84)
in our Lagrangian given in Eq. (6.75), then the equations of motion that we obtain are
163
1 1 1 12 13 1 1 12 13 1
2 2 12 2 23 2 12 2 23 2
3 3 13 23 3 3 13 23 3 3
0 0 2 2 0
0 0 2 2 0
0 0 2 2 0
m x m b b x k b b x
m x b m b x b k b x
m x b b m x b b k x
(6.85)
or equivalently,
1 1 1 12 13 1 1 12 13 1
2 2 12 2 23 2 12 2 23 2
3 3 13 23 3 3 13 23 3 3
0 0 2 2 0
0 0 2 2 0
0 0 2 2 0
m x m b b x k b b x
m x b m b x b k b x
m x b b m x b b k x
(6.86)
in which the damping and stiffness matrices are now both symmetric if we set
13
0 b .
We now generalize this observation to general n-DOF systems. Our aim is to find the
appropriate matrices M , K , and B which when inserted in the Lagrangian given in Eq.
(6.75) will result in the Euler-Lagrange equation given in Eq. (6.79) that is of the form
, Mx Bx Kx 0 (6.87)
where M is a positive definite diagonal matrix, and the matrices B and K are both
symmetric matrices, such as in Eq. (6.86).
Consider the diagonal matrices
1 2
, , ,
n
diag m m m M and
1 2
, , ,
n
diag k k k K
where 2 n . We propose to: (i) change the signs of some of the elements of these
matrices and (ii) provide a procedure to make the matrix M positive definite, and the
matrices B and K symmetric. We do this as follows. If we place negative signs on
i
m ,
j
m ,
k
m , ,
i
k ,
j
k ,
k
k , , , , 1,2, , , and i j k n i j k , then the elements of
the skew-symmetric matrix B which is given by
164
12 13 1
12 23 2
13 23 3
1 2 3
0
0
0
0
n
n
n
n n n
b b b
b b b
b b b
b b b
B
(6.88)
should have its elements altered by the following rule:
(1) The elements are set so that 0
ij
b , 0
ik
b , 0
jk
b , and
(2) After deleting the ith, jth, kth, rows and ith, jth, kth, columns of the B
matrix in Eq. (6.88), the remaining elements of B are set to zero.
Clearly, if we want to place only one negative sign, say, on
i
m and
i
k ( 1,2, , i n ), then
only the second rule (2) applies, since B is skew-symmetric. Also, changing the signs of
all the
i
m ’s and
i
k ’s ( 1,2, , i n ), i.e.,
1 2
, , ,
n
m m m and
1 2
, , ,
n
k k k , will result in
B 0 , a case already considered in Eqs. (6.80) and (6.81), though more generally.
For example, in a 4-DOF system if we have
12 13 14 1 1
12 23 24 2 2
13 23 34 3 3
14 24 34 4 4
0 0 0 0 0 0 0
0 0 0 0 0 0 0
, , ,
0 0 0 0 0 0 0
0 0 0 0 0 0 0
b b b m k
b b b m k
b b b m k
b b b m k
M K B
(6.89)
that is, 2 i and 3 j in this case, then we should choose the elements of the B matrix
by the rule:
(1)
23
0 b , and
(2) After we delete the 2nd and 3rd rows and columns, the remaining elements should
be set to zero, i.e., set
14
0 b .
165
In brief, we should choose the following B matrix:
12 13
12 24
13 34
24 34
0 0
0 0
.
0 0
0 0
b b
b b
b b
b b
B (6.90)
Now using the Lagrangian given by Eq. (6.75), with M and K defined in Eq. (6.89)
and B defined in Eq. (6.90), the equations of motion given by Eq. (6.79) become
1 12 13 1 12 13 1 1 1
12 2 24 12 2 24 2 2 2
13 3 34 13 3 34 3 3 3
24 34 4 24 4 4 4
2 2 0 0 0 0 0
2 0 2 0 0 0 0
2 0 2 0 0 0 0
0 2 2 0 0 0 0
m b b k b b m x x
b m b b k b m x x
b m b b k b m x x
b b m b b m x x
1
2
3
34 4 4
0
0
0
0
x
x
x
k x
(6.91)
or
1 12 13 1 12 13 1 1 1
12 2 24 12 2 24 2 2 2
13 3 34 13 3 34 3 3 3
24 34 4 24 34 4 4 4 4
2 2 0 0 0 0 0
2 0 2 0 0 0 0
2 0 2 0 0 0 0
0 2 2 0 0 0 0
m b b k b b m x x
b m b b k b m x x
b m b b k b m x x
b b m b b k m x x
1
2
3
4
0
0
,
0
0
x
x
x
x
(6.92)
and Eq. (6.92) has symmetric damping and stiffness matrices as well as a positive definite
mass matrix. The corresponding Lagrangian is given by Eq. (6.75).
More generally, when placing negative signs on
i
m ,
j
m ,
k
m , ,
i
k ,
j
k ,
k
k , ,
(and following the procedure described earlier) we have 2
n
different choices, each of
which yields the corresponding equations of motion in the form of Eq. (6.87). However,
by symmetry half of them are equivalent. We can extend this symmetry to the situation
where we consider placing no negative signs on any of the , , 1, 2,. . .,
i i
m k i n (and
following the procedure) to be equivalent to placing negative signs on every
, , 1, 2,. . .,
i i
m k i n , which yields, as mentioned before, B 0 . Due to this symmetry
166
we therefore have
1
2
n
differently ‘structured’ systems of equations of motion of the form
of Eq. (6.87) (each with differently structured matrices B and K ) for which the
corresponding Lagrangians can be obtained by using Eq. (6.75). Furthermore, for each
system structure one can use any (permissible) parameter values. Thus, for large n, one
generates Lagrangians for numerous ordered direct representations of different n-DOF
damped linear systems.
As a special case, in an n-DOF system ( 2 n ) if one chooses
1 1
1
2 2
1 2
3 3
2 3
4 4
3
1 1
0 0 0 0 0 0 0 0
0 0 0 0
0 0 0 0 0 0 0 0
0 0 0
0 0 0 0 0 0 0 0
0 0 0
, ,
0 0 0 0 0 0 0 0
0 0 0 0
0 0 0 0 0 0 0 0 0 1 0 0 0 0 1
n n
n n
m k
b
m k
b b
m k
b b
m k
b
m k
M K B
,
(6.93)
using the matrices given by Eq. (6.93) in the Lagrangian given by Eq. (6.75), the Euler-
Lagrange equation of motion given by Eq. (6.79) yields
1 1 1 1 1
2 2 1 2 2 2
3 3 2 3 3 3
4 4 3 4 4
0 0 0 0 2 0 0 0
0 0 0 0 2 2 0 0
0 0 0 0 0 2 2 0
0 0 0 0 0 0 2 0
0 0 0 0 0 0 0 0
n n n n
m x m b x
m x b m b x
m x b m b x
m x b m x
m x m x
1 1 1
1 2 2 2
2 3 3 3
3 4 4
0 0 0 0
0 0 0
0 0 0
,
0 0 0 0
0 0 0 0 0
n n
k b x
b k b x
b k b x
b k x
k x
(6.94)
which is the equation for an MDOF system with tridiagonal symmetric damping and
stiffness matrices of the form of Eq. (6.87). It should be noted that because the
Lagrangian has negative quantities in the mass matrix, it is not a physical Lagrangian.
Yet, it can give the equations of motion of certain systems consisting of a chain of masses
in which each mass is connected to its neighbors by linear dashpots and springs. The
167
strength of the procedure introduced in this section is that it totally bypasses the
Helmholtz conditions which are near-impossible to solve for arbitrarily large, finite
values of n.
To exemplify what has been discussed so far, we consider a 3-DOF system described
by Eq. (6.86) with
13
0 b , where M , B , and K are given by
1 12 1
2 12 23 2
3 23 3
0 0 0 0 0 0
0 0 , 0 , 0 0 .
0 0 0 0 0 0
m b k
m b b k
m b k
M B K (6.95)
Using the matrices given in Eq. (6.95) in the Lagrangian function given by Eq. (6.75), the
resulting Euler-Lagrange equations of motion are given by Eq. (6.86) which is of the
form of Eq. (6.87), where
1 1 12 1 12
2 12 2 23 12 2 23
3 23 3 23 3
0 0 2 0 0
0 0 , 2 2 , .
0 0 0 2 0
m m b k b
m b m b b k b
m b m b k
M B K
(6.96)
In the theory of linear vibrations, one often considers proportionally damped systems.
One can then particularize the damping matrix B of this system to have proportional
damping so that
, B M K (6.97)
where and are some constants. Substituting Eq. (6.96) into Eq. (6.97), we have the
following relationships:
2 2
, , 1,2,3,
i
i
k
i
m
(6.98)
168
that is, for proportional damping to be possible, each ratio /
i i
k m (see Eq. (6.95)) should
be the same. It is straightforward to verify that Eq. (6.98) still holds for general n-DOF (
1,2, , i n ) tridiagonal systems whose equations of motion are described by Eq. (6.94).
Now, consider a mass-spring-damper system of the type shown in Fig. 6-3 that appears
commonly in the theory of linear vibrations.
Fig. 6-3: General 2-DOF mass-spring-damper system
The equations of motion of the system are
1 1 1
2 2 2
3 3 3
1 0 0 2 1 0 2 1 0 0
0 1 0 1 2 1 1 2 1 0 ,
0 0 1 0 1 2 0 1 2 0
x x x
m x b x k x
x x x
(6.99)
with 0 b . (When 0 b , the forces can be derived from a potential, and a Lagrangian
for the system can be easily found.) We note the tridiagonal form of the damping and
stiffness matrices, and Eq. (6.99) describes a proportionally damped system with 0
and / . b k For this system to be described by the Lagrangian given in Eq. (6.75), in
which the matrices M , B , and K are specified by Eq. (6.95), we require and to
169
also satisfy the conditions given in Eqs. (6.98). Hence, 2 / k b , and noting that
i
m m and 2
i
k k ( 1,2,3 i ) we require
.
k b
b m
(6.100)
Application of this result to n-degree-of-freedom tridiagonal systems follows easily
by extension. In fact, the condition Eq. (6.100) along with 2 / k b should be also met
for general n-DOF tridiagonal systems of the type given in Eq. (6.99) to be proportionally
damped.
To illustrate this result, we apply it to a proportionally damped system described by
Eq. (6.99) with the parameters 0.64 kg m , 0.8 N s/m b , and 1 N/m k . Noting that
Eq. (6.100) holds, we find
1
2
2.5 s ,
k
b
1
0 s ,
0.8 s,
b
k
1 2 3
0.64 kg, m m m m (6.101a)
12
0.4 N s/m,
2
b
b
23
0.4 N s/m,
2
b
b
1 2 3
2 2 N/m, k k k k (6.101b)
by comparing Eq. (6.99) with Eqs. (6.87) and (6.96). The Lagrangian for this dissipative
3-DOF system given by
2.5 2 2 2 2 2 2
1 2 3 1 2 3 1 2 2 3 1 2 2 3
1 1
2 2
0.32 0.4
t
t
L e
e x x x x x x x x x x x x x x
T T T
x Mx x Kx x Bx
(6.102)
then yields the equations of motion given in Eq. (6.99) for the specific values of m, k, and
b, chosen in this example.
170
As the last application of the results obtained so far, we propose a Jacobi integral that
is conserved at all times for the types of linear MDOF systems considered here (see Eq.
(6.79)). When the Lagrangian does not contain time explicitly (and the actual velocity is a
virtual velocity), i.e., when , , / 0 L t t q q , the Jacobi integral I given by (Pars [54])
:
L
I L
T
q
q
(6.103)
is conserved. The Lagrangian in Eq. (6.75), however, contains time explicitly, but by
using the transformation
2
t
e
y x (6.104)
it becomes
1 1
.
2 2 2 2 2
L
T
T T
y y M y y y Ky y B y y (6.105)
Hence, the Jacobi integral is given by Eq. (6.103) as
1 1
,
2 2 2 2 2
I
T
T T T T
y y M y y y B B y y Ky y By (6.106)
and since the matrix B is assumed to be skew-symmetric in this paper, Eq. (6.106)
simplifies to
1 1
.
2 2 2 2
I
T
T
y y M y y y Ky (6.107)
Rewriting this Jacobi integral in terms of x and x by using Eq. (6.104), we obtain the
conservation law
171
1 1
, , constant.
2 2
t
I t e
T
T
x x x x Mx x Kx (6.108)
Thus, for an MDOF linear system whose equation of motion is given by Eq. (6.79) it is
guaranteed that the function , , I t x x in Eq. (6.108) is conserved at all times. For
example, corresponding to the Lagrangian given in Eq. (6.102) for the dissipative,
proportionally damped 3-DOF system shown in Fig. 6-3 with the parameters given in
Eqs. (6.101), we obtain the conservation law
2.5 2 2 2 2 2 2
1 2 3 1 1 2 2 3 3 1 2 3
, , 0.32 0.8 constant.
t
I t e x x x x x x x x x x x x
x x (6.109)
172
Chapter 7
GENERAL CONCLUSIONS
As for the first topic, starting with an unperturbed, circular reference of the leader
satellite, more and more complicated formation-keeping problems are handled, showing
the progress of my research. First, the problem was solved wholly in the Hill frame and
the case when the reference orbit of the leader is arbitrary is considered. Next, incorrect
initial conditions were handled and the effect of the nonuniform gravity of the
nonspherical Earth was included. Then, it was shown that the explicit solution to the
formation-keeping even with the attitude constraints can be readily obtained, leading to
the separation principle which tells under what conditions orbital control is independent
of the attitude of the satellite. Finally, an additional additive controller based on a
generalization of the concept of sliding surface control was designed to compensate for
model uncertainties.
Regarding the second topic, the existing results about the inverse problem of the
calculus of variations for SDOF systems were extended, and a more systematic way was
derived to get the candidates for nonpotential forces that have their corresponding
integrand functions in Hamilton’s principle. Then, new results for damped, linear MDOF
systems are obtained using extensions of results for SDOF and 2-DOF systems.
173
Conservation laws for these damped MDOF systems are found using the Lagrangians
obtained, and several examples are provided.
The main contributions of this dissertation are the following:
1. We obtain the exact closed-form solution (for the dynamic model assumed) to
the full problem of formation-keeping with both orbital and attitude requirements. Unlike
previous research, we directly start from the nonlinear equation and the exact control
force is easily obtained that captures all the nonlinearities of the nonlinear, non-
autonomous system of differential equations.
2. The newly obtained solution is valid even when the leader satellite is in an
arbitrary orbit and there are initial (insertion) orbital and/or attitude errors. Since the
control force to be applied to the follower satellite(s) is explicitly obtained in closed form
and the method is not computationally intensive, it can be easily used for on-orbit, real-
time control of maneuvers for which the underlying dynamics is highly nonlinear.
3. Since the orbital constraints and the attitude target pointing constraints are
independent of each other (as shown in Appendix A), it is possible to obtain the explicit
form of the solutions, which enables further insights into the dynamics of a multiple-
satellite system. Specifically, it is shown under what conditions the attitude dynamics can
be ignored, or separated away, when designing a controller for exact orbital control,
leading to the statement of the Separation Principle. As a special case, if there are three
174
independent orbital constraints given, the separation principle is always satisfied
throughout the maneuver.
4. It is pointed out why the single constraint (instead of two constraints) for
attitude control of target-pointing fails. Generally speaking, if there are two constraints
1
0 and
2
0 , then we should differentiate each constraint with respect to time to
get constraint equations of the form Eq. (2.6). The two constraints are mathematically
equivalent to
2 2
1 2
0 if they have real values, but this single equation cannot be used
as the constraint equation because differentiation of this single equation only yields
0 q 0 which tells nothing about the constraints.
5. When considering model uncertainties, the generalized control force from the
nominal controller given by the fundamental equation is augmented with an additive
controller based on the generalization of the concept of the sliding surface control. Given
bounds on the uncertainties, the additive controller developed is continuous. The additive
control guarantees that the uncertain system, with unknown parameters, tracks of the
nominal system trajectories within specified error bounds.
6. New integrand functions for some general nonpotential forces are obtained for
SDOF systems and a simple and systematic approach for getting them is developed. It
was shown that the main difficulties lie in performing the necessary integrations
explicitly.
175
7. For a 2-DOF linear system with linear damping, the conditions for the existence
of a Lagrangian are explicitly obtained by solving the Helmholtz conditions. Three
general cases when such Lagrangians are guaranteed to exist are obtained, depending on
the parameter values of the coupled linear systems.
8. By using and generalizing results for SDOF systems, a simple procedure that
does not require the use of the Helmholtz conditions and that is easily extended to n-DOF
linear systems, is developed. We specifically include systems that commonly arise in the
theory of linear vibrations--systems with positive definite mass matrices, and symmetric
stiffness and damping matrices. The method yields several new Lagrangians for linear
MDOF systems.
9. Conservation laws for dissipative MDOF systems are also obtained by finding
the corresponding Jacobi integrals.
176
Chapter 8
FUTURE WORK
Although the approaches adopted in this dissertation could help to solve the formation-
keeping problem and the inverse problem of the calculus of variations in some aspects,
there are still many difficulties to surmount and to be improved. As for the first topic, it is
necessary to include other perturbations on both orbit and attitude in space into the
system such as air drag, solar radiation pressure, gravity gradient torque, and so on. The
fact that modelings about these perturbations are not perfect is one more reason an
additional additive controller is necessary. In this dissertation, two controllers were used
and combined that employ the applications of the fundamental equation and the sliding
surface control. Another combination with other existing control techniques such as
feedback linearization, adaptive control, and so on, can be possible, and the comparison
with the current methodology developed in this dissertation would be also be interesting.
Also, the control proposed herein is an exact control, which is sometime very difficult to
achieve in real life. Also, it may necessitate more energy than other control methods that
allow some larger error bounds. These disadvantages of the current control methodology
need to be improved in future work.
For the second topic, the Bolza’s approach was used to get a Lagrangian for a SDOF
system. The main difficulty of this approach lies in the integration given by Eq. (5.9).
177
These integrations are not possible for many nonlinear differential equations whose
Lagrangians have not yet been known (including van der Pol oscillator). One may find a
transformation that makes the integration possible. Also, a combination with other
methods, including the Lie symmetry method, is expected to help solve the problem. For
MDOF systems, we have dealt with only linear systems in this dissertation. A mass-
spring-damper system with nonlinear (ex. quadratic) damping is feasible and expected to
be solved by a similar extension from a SDOF case. As shown in this dissertation, the
inverse problem of the calculus of variations is highly related with finding the first
integrals or the conservations laws of the system. If a Lagrangian is luckily obtained but
the Lagrangian constains time explicitly, then we cannot use the Jacobi integrals. In this
case, we need to transform the equations to the ones that are directly integrable or that
have Lagrangians which does not contain time explicitly. More generally, a Killing’s
equation, derived from the symmetry of the system by Noether’s theorem, can be used to
find conservation laws and this is still an open and interesting problem.
178
Bibliography
[1] Hill, G. W., “Researches in the Lunar Theory”, American Journal of Mathematics,
Vol. 1, No. 1, 1878, pp. 5-26.
[2] Clohessy, W. H. and Wiltshire, R. S., “Terminal Guidance System for Satellite
Rendezvous,” Journal of Aerospace Sciences, Vol. 27, No. 8, 1960, pp. 653-658, 674.
[3] Tschauner, J. and Hempel, P., “Rendezvous zu Einemin Elliptischer Bahn
Umlaufenden Ziel,” Astronautica Acta, Vol. 11, No. 2, 1965, pp. 104-109.
[4] Schweighart, S. A. and Sedwick, R. J., “High-Fidelity Linearized J2 Model for
Satellite Formation Flight,” Journal of Guidance, Control, and Dynamics, Vol. 25, No. 6,
2002, pp. 1073-1080.
[5] Yan, Q., Kapila, V., and Sparks, A. G., “Pulse-Based Periodic Control for Spacecraft
Formation Flying,” American Control Conference, Chicago, Illinois, June 2000, pp. 374-
378.
[6] Sparks, A., “Satellite Formationkeeping Control in the Presence of Gravity
Perturbations,” American Control Conference, Chicago, Illinois, June 2000, pp. 844-848.
[7] Vaddi, S. S. and Vadali, S. R., “Linear and Nonlinear Control Laws for Formation
Flying,” The 13th AAS/AIAA Space Flight Mechanics Meeting, Paper AAS 03-109,
Ponce, Puerto Rico, February 2003.
[8] Milam, M. B., Petit, N., and Murray, R. M., “Constrained Trajectory Generation for
Micro-Satellite Formation Flying,” AIAA Guidance, Navigation, and Control Conference
and Exhibit, Paper AIAA 2001-4030, Montreal, Canada, August 2001.
[9] De Queiroz, M. S., Kapila, V., and Yan, Q., “Adaptive Nonlinear Control of Multiple
Spacecraft Formation Flying,” Journal of Guidance, Control, and Dynamics, Vol. 23,
No. 3, 2000, pp. 385-390.
[10] Kristiansen, R., Grotli, E. I., Nicklasson, P. J., and Gravdahl, J. T., “A Model of
Relative Position and Attitude in a Leader-Follower Spacecraft Formation,”
Scandinavian Conference on Simulation and Modeling, Trondheim, Norway, June 2005.
[11] Long, M. R. and Hall, C. D., “Attitude Tracking Control for Spacecraft Formation
Flying,” Flight Mechanics Symposium, Goddard Space Flight Center, May 1999.
179
[12] Ahn, C. and Kim, Y., “Point Targeting of Multisatellites via a Virtual Structure
Formation Flight Scheme,” Journal of Guidance, Control, and Dynamics, Vol. 32, No. 4,
2009, pp. 1330-1344.
[13] Segal, S. and Gurfil, P., “Effect of Kinematic Rotation-Translation Coupling on
Relative Spacecraft Translational Dynamics,” Journal of Guidance, Control, and
Dynamics, Vol. 32, No. 3, 2009, pp. 1045-1050.
[14] Udwadia, F. E. and Kalaba, R. E., “What Is the General Form of the Explicit
Equations of Motion for Constrained Mechanical Systems,” Journal of Applied
Mechanics, Vol. 69, No. 3, 2002, pp. 335-339.
[15] Udwadia, F. E., “Equations of Motion for Mechanical Systems: A Unified
Approach,” International Journal of Non-Linear Mechanics, Vol. 31, No. 6, 1996, pp.
951-958.
[16] Udwadia, F. E., “Nonideal Constraints and Lagrangian Dynamics,” Journal of
Aerospace Engineering, Vol. 13, No. 1, 2000, pp. 17-22.
[17] Kalaba, R. E. and Udwadia, F. E., “Equations of Motion for Nonholonomic,
Constrained Dynamical Systems via Gauss’s Principle,” Journal of Applied Mechanics,
Vol. 60, No. 3, 1993, pp. 662-668.
[18] Udwadia, F. E., “A New Perspective on the Tracking Control of Nonlinear
Structural and Mechanical Systems,” Proceedings of the Royal Society A, Vol. 459, 2003,
pp. 1783-1800.
[19] Udwadia, F. E., “Equations of Motion for Constrained Multibody Systems and Their
Control,” Journal of Optimization Theory and Applications, Vol. 127, No. 3, 2005, pp.
627-638.
[20] Udwadia, F. E., “Optimal Tracking Control of Nonlinear Dynamical Systems,”
Proceedings of the Royal Society A, Vol. 464, 2008, pp. 2341-2363.
[21] Udwadia, F. E. and Kalaba, R. E., “On Motion,” Journal of the Franklin Institute,
Vol. 330 , No. 3, 1993, pp. 571-577.
[22] Udwadia, F. E. and Kalaba, R. E., “A New Perspective on Constrained Motion,”
Proceedings of the Royal Society A, Vol. 439, 1992, pp. 407-410.
180
[23] Udwadia, F. E. and Kalaba, R. E., “An Alternative Proof of the Greville Formula,”
Journal of Optimization Theory and Application, Vol. 94, No. 1, 1997, pp. 23-28.
[24] Schutte, A. D. and Dooley, B. A., “Constrained Motion of Tethered Satellites,”
Journal of Aerospace Engineering, Vol. 18, No. 4, 2005, pp. 242-250.
[25] Alizadeh, I. and Villac, B., “Applications of Constraint Stabilization to Low-Thrust
Mission Design,” The 19th AAS/AIAA Space Flight Mechanics Meeting, Paper AAS 09-
146, Savannah, Georgia, February 2009.
[26] Udwadia, F. E., Schutte, A. D., and Lam, T., “Formation Flight of Multiple Rigid
Body Spacecraft,” The 48th AIAA/ASME/ASCE/AHS/ASC Structures, Structural
Dynamics, and Materials Conference, Paper AIAA 2007-2391, Honolulu, Hawaii, April
2007.
[27] De Queirox, M. S., Yan, Q., Yang, G., and Kaplia, V., “Global Output Feedback
Tracking Control of Spacecraft Formation Flying with Parametric Uncertainty,”
Proceedings of the 38
th
IEEE Conference on Decision and Control, 1999, Vol. 1, pp.
584-589.
[28] Wang, Z., Khorrami, F., and Grossman, W., “Robust Adaptive Control of
Formationkeeping for a Pair of Satellites,” Proceedings of the 2000 American Control
Conference, 2000, Vol. 2, pp. 834-838.
[29] Gautam, A., Soh, Y. C., and Chu, Y. –C., “Robust Predictive Control of Satellite
Formations,” Proceedings of the 10
th
International Conference on Control, Automation,
Robotics and Vision, 2008, pp. 844-849.
[30] Bolza, O., Vorlesungen über Varationsrechnung, Koehler und Amelang, Leipzig,
Germany, 1931.
[31] Douglas, J., “Theorems in the Inverse Problem in the Calculus of Variations,”
Proceedings of the National Academy of Sciences, Vol. 26, 1940, pp. 215–221.
[32] Douglas, J., “Solution of the Inverse Problem of the Calculus of Variations,”
Transactions of the American Mathematical Society, Vol. 50, 1941, pp. 71-128.
[33] Leitmann, G., “Some Remarks on Hamilton’s Principle,” Journal of Applied
Mechanics, Vol. 30, 1963, pp. 623-625.
181
[34] He, J. -H., “Variational Principles for Some Nonlinear Partial Differential Equations
with Variable Coefficients,” Chaos, Solitons, and Fractals, Vol. 19, 2004, pp. 847-851.
[35] Darboux, G., Leçons sur la Théorie Générale des Surfaces, Vol. 3, Gauthier-Villars,
Paris, 1894.
[36] Helmholtz, H. v., “Uber die physikalische Bedeutung des Princips der kleinsten
Wirkung,” Journal für die reine angewandte Mathematik, Vol. 100, 1887, pp. 137-166.
[37] Santilli, R. M., Foundations of Theoretical Mechanics I: The Inverse Problem in
Newtonian Mechanics, Springer-Verlag, New York, 1978.
[38] Hojman, S. and Ramos, S., “Two-Dimensional s-Equivalent Lagrangians and
Separability,” Journal of Physics A: Mathematical and General, Vol. 15, 1982, pp. 3475-
3480.
[39] Mestdag, T., Sarlet, W., and Crampin, M., “The Inverse Problem for Lagrangian
Systems with Certain Non-conservative Forces,” Differential Geometry and Its
Applications, Vol. 29, 2011, pp. 55-72.
[40] Vallado, D. A., Fundamentals of Astrodynamics and Applications, 2nd ed., The
Space Technology Library, El Segundo, 2002.
[41] Baumgarte, J., “Stabilization of Constraints and Integrals of Motion in Dynamical
Systems,” Computer Methods in Applied Mechanics and Engineering, Vol. 1, No. 1,
1972, pp. 1-16.
[42] Udwadia, F. E. and Schutte, A. D., “An Alternative Derivation of the Quaternion
Equations of Motion for Rigid-Body Rotational Dynamics,” Journal of Applied
Mechanics, Vol. 77, No. 4, 2010, pp. 044505-1 - 044505-4.
[43] Sabol, C., Burns, R., and McLaughlin, C. A., “Satellite Formation Flying Design and
Evolution,” Journal of Spacecraft and Rockets, Vol. 38, No. 2, 2001, pp. 270-278.
[44] Udwadia, F. E. and Kalaba, R. E., Analytical Dynamics: A New Approach,
Cambridge University Press, New York, 2008.
[45] Busse, F., Inalhan, G., and How, J. P., “Project ORION: Carrier Phase Differential
GPS Navigation for Formation Flying,” Annual AAS Rocky Mountain Conference,
Breckenridge, Colorado, February 2000, pp. 197-212.
182
[46] How, J. P. and Tillerson M., “Analysis of the Impact of Sensor Noise on Formation
Flying Control,” American Control Conference, Arlington, Virginia, June 2001, pp.
3986-3991.
[47] Henderson, H. V. and Searle, S. R., “On Deriving the Inverse of a Sum of Matrices,”
SIAM Review, Vol. 23, No. 1, 1981, pp. 53-60.
[48] Utkin, V. I., “Variable Structure with Sliding Mode-A Survey,” IEEE Transactions
of Automatic Control, Vol. 22, No. 2, 1977, 212-222.
[49] Tsypkin and Ya, Z., “Teoriya Releinykh System Avtomati Cheskogo Regulirovaniya
(Theory of Switching Control Systems),” Gostekhizdat, Moscow, 1995.
[50] Edwards, C. and Spurgeon, S., “Sliding Mode Control: Theory and Applications,”
Taylor and Francis, London, 1999
[51] Khalil, H. K., Nonlinear Systems, Prentice-Hall, Upper Saddle River, New Jersey,
2002.
[52] Leitmann, G., The Calculus of Variations and Optimal Control: An Introduction,
Springer, 1986, New York.
[53] Musielak, Z. E., 2008, “Standard and Non-Standard Lagrangians for Dissipative
Dynamical Systems With Variable Coefficients,” Journal of Physics A: Mathematical
Theories, Vol. 41, p. 055205.
[54] Ray, J. R., “Lagrangians and Systems They Describe - How Not to Treat Dissipation
in Quantum Mechanics,” American Journal of Physics, Vol. 47, 1979, pp. 626-629.
[55] Pars, L. A., A Treatise on Analytical Dynamics, Ox Bow Press, 1972.
[56] Horn, R. A. and Johnson, C. R., Matrix Analysis, Cambridge University Press, New
York, 1985.
183
Appendix A
In this appendix we give explicit expressions for the elements of the matrices R , R
, and
R
. Let us first define the following expressions that will be used frequently:
1/2
2 2 2
L L L L
r X Y Z ,
L L L L L L
L
L
X X Y Y Z Z
r
r
,
2 2 2
2
1 L L L L L L L L L L
L
L
L L L L L L L
r X Y Z X X Y Y Z Z
r
r
r X X Y Y Z Z
,
1/2
2 2 2
L L L L L L L L L L L L L
h Y Z Z Y Z X X Z X Y Y X
,
1
L L L L L L L L L L L L L L L L L L L L L L L L L
L
h Y Z Z Y Y Z Z Y Z X X Z Z X X Z X Y Y X X Y Y X
h
,
and,
2
2
2
2
1
L L L L L L L L L L L L L L L L
L L L L L L L L L L L L L L L L L
L
L L L L L L L L L L L L L L L L
L
Y Z Z Y Y Z Z Y Y Z Y Z Z Y Z Y
h Z X X Z Z X X Z Z X Z X X Z X Z
h
X Y Y X X Y Y X X Y X Y Y X Y X
h
L L L L L L L L L L L L L L L L
L
L L L L L L L L
Y Z Z Y Y Z Z Y Z X X Z Z X X Z
h
X Y Y X X Y Y X
.
Then, the rotation matrix
T
R S , and its derivatives R
, R
are of the form
11 12 13
21 22 23
31 32 33
R R R
R R R
R R R
R
,
11 12 13
21 22 23
31 32 33
R R R
R R R
R R R
R
,
11 12 13
21 22 23
31 32 33
R R R
R R R
R R R
R
,
where each element of these matrices is given as follows:
11
L
L
X
R
r
,
11 2
L L L L
L
r X X r
R
r
,
11 3
1
2
L L L L L L L L L L
L
R r r X X r r r X X r
r
,
12
L
L
Y
R
r
,
12 2
L L L L
L
r Y Y r
R
r
,
12 3
1
2
L L L L L L L L L L
L
R r r Y Y r r r Y Y r
r
,
184
13
L
L
Z
R
r
,
13 2
L L L L
L
r Z Z r
R
r
,
13 3
1
2
L L L L L L L L L L
L
R r r Z Z r r r Z Z r
r
,
21
1
L L L L L L L L L L
L L
R Z Z X X Z Y X Y Y X
h r
,
21 2 2
1
L L L L L L L L L L L L L L L L L L L L L L
L L
L L L L L L L L L L L L L L
h r Z Z X X Z Z Z X X Z Y X Y Y X Y X Y Y X
R
h r
h r h r Z Z X X Z Y X Y Y X
,
21 3 3
2
2
1
L L L L L L L L L L L L L L L L L L L L L L L L
L L L L L L L L L L L L L L L L L L L
L L
L L L L L L L L L L L L L
L L
h r h r Z Z X X Z Z Z X X Z Y X Y Y X Y X Y Y X
Z Z X X Z Z Z X X Z Z Z X Z X X Z X Z
h r
Y X Y Y X Y X Y Y X Y X h r
R
h r
2
L L L L L L L L
L L L L L L L L L L L L L L L L
L L L L L L L L L L L L L L L L L L L L L L L L
Y X Y Y X Y X
h r h r h r Z Z X X Z Y X Y Y X
h r h r Z Z X X Z Z Z X X Z Y X Y Y X Y X Y Y X
2
L L L L L L L L L L L L L L L L L L L L L L
L L L L
L L L L L L L L L L L L L L
h r Z Z X X Z Z Z X X Z Y X Y Y X Y X Y Y X
h r h r
h r h r Z Z X X Z Y X Y Y X
,
22
1
L L L L L L L L L L
L L
R X X Y Y X Z Y Z Z Y
h r
,
22 2 2
1
L L L L L L L L L L L L L L L L L L L L L L
L L
L L L L L L L L L L L L L L
h r X X Y Y X X X Y Y X Z Y Z Z Y Z Y Z Z Y
R
h r
h r h r X X Y Y X Z Y Z Z Y
,
22 3 3
2
2
1
L L L L L L L L L L L L L L L L L L L L L L L L
L L L L L L L L L L L L L L L L L L L
L L
L L L L L L L L L L L L L
L L
h r h r X X Y Y X X X Y Y X Z Y Z Z Y Z Y Z Z Y
X X Y Y X X X Y Y X X X Y X Y Y X Y X
h r
Z Y Z Z Y Z Y Z Z Y Z Y h r
R
h r
2
L L L L L L L L
L L L L L L L L L L L L L L L L
L L L L L L L L L L L L L L L L L L L L L L L L
Z Y Z Z Y Z Y
h r h r h r X X Y Y X Z Y Z Z Y
h r h r X X Y Y X X X Y Y X Z Y Z Z Y Z Y Z Z Y
2
L L L L L L L L L L L L L L L L L L L L L L
L L L L
L L L L L L L L L L L L L L
h r X X Y Y X X X Y Y X Z Y Z Z Y Z Y Z Z Y
h r h r
h r h r X X Y Y X Z Y Z Z Y
,
23
1
L L L L L L L L L L
L L
R Y Y Z Z Y X Z X X Z
h r
,
23 2 2
1
L L L L L L L L L L L L L L L L L L L L L L
L L
L L L L L L L L L L L L L L
h r Y Y Z Z Y Y Y Z Z Y X Z X X Z X Z X X Z
R
h r
h r h r Y Y Z Z Y X Z X X Z
,
185
23 3 3
2
2
1
L L L L L L L L L L L L L L L L L L L L L L L L
L L L L L L L L L L L L L L L L L L L
L L
L L L L L L L L L L L L L
L L
h r h r Y Y Z Z Y Y Y Z Z Y X Z X X Z X Z X X Z
Y Y Z Z Y Y Y Z Z Y Y Y Z Y Z Z Y Z Y
h r
X Z X X Z X Z X X Z X Z h r
R
h r
2
L L L L L L L L
L L L L L L L L L L L L L L L L
L L L L L L L L L L L L L L L L L L L L L L L L
X Z X X Z X Z
h r h r h r Y Y Z Z Y X Z X X Z
h r h r Y Y Z Z Y Y Y Z Z Y X Z X X Z X Z X X Z
2
L L L L L L L L L L L L L L L L L L L L L L
L L L L
L L L L L L L L L L L L L L
h r Y Y Z Z Y Y Y Z Z Y X Z X X Z X Z X X Z
h r h r
h r h r Y Y Z Z Y X Z X X Z
,
31
L L L L
L
Y Z Z Y
R
h
,
31 2
1
L L L L L L L L L L
L
R h Y Z Z Y h Y Z Z Y
h
,
31 3
1
2
L L L L L L L L L L L L L L L
L
L L L L L L L L L L L
h h Y Z Y Z Z Y Z Y h Y Z Z Y
R
h
h h Y Z Z Y h Y Z Z Y
,
32
L L L L
L
Z X X Z
R
h
,
32 2
1
L L L L L L L L L L
L
R h Z X X Z h Z X X Z
h
,
32 3
1
2
L L L L L L L L L L L L L L L
L
L L L L L L L L L L L
h h Z X Z X X Z X Z h Z X X Z
R
h
h h Z X X Z h Z X X Z
,
33
L L L L
L
X Y Y X
R
h
,
33 2
1
L L L L L L L L L L
L
R h X Y Y X h X Y Y X
h
,
and
33 3
1
2
L L L L L L L L L L L L L L L
L
L L L L L L L L L L L
h h X Y X Y Y X Y X h X Y Y X
R
h
h h X Y Y X h X Y Y X
.
186
Appendix B
This appendix provides a Maple code that generates the time derivatives of the
acceleration vectors of a satellite in ITRF under the nonuniform gravity perturbation of
the nonspherical Earth. After running this code, one can directly use the result in Matlab.
> restart; with(orthopoly);
> LAssocPoly[0, 0] := 1; LAssocPoly[0, 1] := 0; LAssocPoly[2, 0] := 1/2*(3*z(t)^2/r^2-
1); LAssocPoly[2, 1] := 3*z(t)*sqrt(1-z(t)^2/r^2)/r; LAssocPoly[2, 2] := 3*(1-z(t)^2/r^2);
LAssocPoly[2, 3] := 0; LAssocPoly[3, 0] := (5/2)*z(t)^3/r^3-(3/2)*z(t)/r; LAssocPoly[3,
1] := (1/2*(15*z(t)^2/r^2-3))*sqrt(1-z(t)^2/r^2); LAssocPoly[3, 2] := (15*(1-
z(t)^2/r^2))*z(t)/r; LAssocPoly[3, 3] := 15*(1-z(t)^2/r^2)^(3/2); LAssocPoly[3, 4] := 0;
LAssocPoly[4, 0] := 1/8*(35*z(t)^4/r^4-30*z(t)^2/r^2+3); LAssocPoly[4, 1] :=
(5/2)*sqrt(1-z(t)^2/r^2)*(7*z(t)^3/r^3-3*z(t)/r); LAssocPoly[4, 2] := (15/2*(1-
z(t)^2/r^2))*(7*z(t)^2/r^2-1); LAssocPoly[4, 3] := 105*z(t)*(1-z(t)^2/r^2)^(3/2)/r;
LAssocPoly[4, 4] := 105*(1-z(t)^2/r^2)^2; LAssocPoly[4, 5] := 0; lmax := 5; for ll from
0 to lmax do for mm from 0 to ll do lll := ll; mmm := mm; nn := lll+mmm; PP := (1-
v^2)^((1/2)*mmm)*(diff((v^2-1)^lll, [`$`(v, nn)]))/(2^lll*factorial(lll)); ff := subs(v =
sin(phi), PP); ff1 := subs(csgn(cos(phi)) = 1, simplify(ff)); LAssocPoly[lll, mm] :=
subs({cos(phi) = sqrt(1-z(t)^2/r^2), sin(phi) = z(t)/r}, ff1) end do; LAssocPoly[ll, ll+1] :=
0 end do;
print(`output redirected...`); # input placeholder
> deltaUbydeltar1 := -mu*x(t)*(sum(sum((R/r)^l*(l+1)*LAssocPoly[l, m]*(C(l,
m+1)*cos(m*MTM[atan](y(t)/x(t)))+S(l, m+1)*sin(m*MTM[atan](y(t)/x(t)))), m = 0 ..
l), l = 2 .. lmax))/r^3;
print(`output redirected...`); # input placeholder
> deltaUbydeltar11 := subs(r = sqrt(x(t)^2+y(t)^2+z(t)^2), deltaUbydeltar1); X1 :=
simplify(diff(deltaUbydeltar11, t), short);
print(`output redirected...`); # input placeholder
print(`output redirected...`); # input placeholder
> deltaUbydeltaphi1 := mu*z(t)*x(t)*(sum(sum((R/r)^l*(LAssocPoly[l, m+1]-
m*tan(MTM[asin](z(t)/r))*LAssocPoly[l, m])*(C(l,
187
m+1)*cos(m*MTM[atan](y(t)/x(t)))+S(l, m+1)*sin(m*MTM[atan](y(t)/x(t)))), m = 0 ..
l), l = 2 .. lmax))/(r^3*sqrt(x(t)^2+y(t)^2));
print(`output redirected...`); # input placeholder
> deltaUbydeltaphi11 := subs(r = sqrt(x(t)^2+y(t)^2+z(t)^2), deltaUbydeltaphi1); X2 :=
simplify(-(diff(deltaUbydeltaphi11, t)), short);
print(`output redirected...`); # input placeholder
> deltaUbydeltalambda1 := mu*y(t)*(sum(sum((R/r)^l*m*LAssocPoly[l, m]*(S(l,
m+1)*cos(m*MTM[atan](y(t)/x(t)))-C(l, m+1)*sin(m*MTM[atan](y(t)/x(t)))), m = 0 .. l),
l = 2 .. lmax))/(r*(x(t)^2+y(t)^2));
print(`output redirected...`); # input placeholder
> deltaUbydeltalambda11 := subs(r = sqrt(x(t)^2+y(t)^2+z(t)^2), deltaUbydeltalambda1);
print(`output redirected...`); # input placeholder
> X3 := simplify(diff(-deltaUbydeltalambda11, t), short);
print(`output redirected...`); # input placeholder
> removeTime := [diff(x(t), t) = x, diff(y(t), t) = y, diff(z(t), t) = z, x(t) = x, y(t) = y, z(t) =
z];
print(`output redirected...`); # input placeholder
> XX1 := subs(removeTime, X1);
print(`output redirected...`); # input placeholder
> XX2 := subs(removeTime, X2);
> XX3 := subs(removeTime, X3);
> XXX := simplify(XX1+XX2+XX3, short);
print(`output redirected...`); # input placeholder
> deltaUbydeltar2 := -mu*y(t)*(sum(sum((R/r)^l*(l+1)*LAssocPoly[l, m]*(C(l,
m+1)*cos(m*MTM[atan](y(t)/x(t)))+S(l, m+1)*sin(m*MTM[atan](y(t)/x(t)))), m = 0 ..
l), l = 2 .. lmax))/r^3;
%;
> deltaUbydeltar22 := subs(r = sqrt(x(t)^2+y(t)^2+z(t)^2), deltaUbydeltar2);
> Y1 := simplify(diff(deltaUbydeltar22, t), short);
print(`output redirected...`); # input placeholder
> deltaUbydeltaphi2 := mu*z(t)*y(t)*(sum(sum((R/r)^l*(LAssocPoly[l, m+1]-
m*tan(MTM[asin](z(t)/r))*LAssocPoly[l, m])*(C(l,
m+1)*cos(m*MTM[atan](y(t)/x(t)))+S(l, m+1)*sin(m*MTM[atan](y(t)/x(t)))), m = 0 ..
l), l = 2 .. lmax))/(r^3*sqrt(x(t)^2+y(t)^2));
> deltaUbydeltaphi22 := subs(r = sqrt(x(t)^2+y(t)^2+z(t)^2), deltaUbydeltaphi2);
> Y2 := simplify(-(diff(deltaUbydeltaphi22, t)), short);
print(`output redirected...`); # input placeholder
188
> deltaUbydeltalambda2 := mu*x(t)*(sum(sum((R/r)^l*m*LAssocPoly[l, m]*(S(l,
m+1)*cos(m*MTM[atan](y(t)/x(t)))-C(l, m+1)*sin(m*MTM[atan](y(t)/x(t)))), m = 0 .. l),
l = 2 .. lmax))/(r*(x(t)^2+y(t)^2));
> deltaUbydeltalambda22 := subs(r = sqrt(x(t)^2+y(t)^2+z(t)^2), deltaUbydeltalambda2);
> Y3 := simplify(diff(deltaUbydeltalambda22, t), short);
print(`output redirected...`); # input placeholder
> YY1 := subs(removeTime, Y1);
print(`output redirected...`); # input placeholder
> YY2 := subs(removeTime, Y2);
> YY3 := subs(removeTime, Y3);
> YYY := simplify(YY1+YY2+YY3, short);
print(`output redirected...`); # input placeholder
> deltaUbydeltar3 := -mu*z(t)*(sum(sum((R/r)^l*(l+1)*LAssocPoly[l, m]*(C(l,
m+1)*cos(m*MTM[atan](y(t)/x(t)))+S(l, m+1)*sin(m*MTM[atan](y(t)/x(t)))), m = 0 ..
l), l = 2 .. lmax))/r^3; deltaUbydeltar33 := subs(r = sqrt(x(t)^2+y(t)^2+z(t)^2),
deltaUbydeltar3);
> Z1 := simplify(diff(deltaUbydeltar33, t), short);
print(`output redirected...`); # input placeholder
> deltaUbydeltaphi3 := mu*sqrt(x(t)^2+y(t)^2)*(sum(sum((R/r)^l*(LAssocPoly[l, m+1]-
m*tan(MTM[asin](z(t)/r))*LAssocPoly[l, m])*(C(l,
m+1)*cos(m*MTM[atan](y(t)/x(t)))+S(l, m+1)*sin(m*MTM[atan](y(t)/x(t)))), m = 0 ..
l), l = 2 .. lmax))/r^3;
%;
> deltaUbydeltaphi33 := subs(r = sqrt(x(t)^2+y(t)^2+z(t)^2), deltaUbydeltaphi3);
> Z2 := simplify(diff(deltaUbydeltaphi33, t), short);
print(`output redirected...`); # input placeholder
> ZZ1 := subs(removeTime, Z1);
print(`output redirected...`); # input placeholder
> ZZ2 := subs(removeTime, Z2);
print(`output redirected...`); # input placeholder
>
> ZZZ := simplify(ZZ1+ZZ2, short);
print(`output redirected...`); # input placeholder
> with(CodeGeneration); Matlab(XXX, resultname = "accX"); Matlab(YYY, resultname
= "accY"); Matlab(ZZZ, resultname = "accZ");
189
Appendix C
Here it is shown that the matrix A in Eq. (3.36) has always full rank as long as
11
A is
full rank. The assumption that
11
A has full rank, is equivalent to assuming that the orbital
tracking requirements for the ith follower satellite are independent of one another.
Let us first show that the matrix
21
A has always full rank, that is, rank of two. From
the definition of
21
A in Eq. (3.33),
.
2
21
1
T
A
T
(C-1)
Since the transformation matrix
1
2
3
T
T T
T
is orthogonal, it is obvious that the two rows of
21
A are linearly independent, which proves that
21
A has full rank.
Next it is shown that the 3 by 4 matrix
22
A has always full rank, namely rank of 3.
Because the matrix E (see Eq. (2.64)) is orthogonal,
rank rank
T
22 22
A A E , and
from Eqs. (3.34) and (2.64), one has
,
2 2
T
T T 22 22 1 22
22 1
T
1 3
A u A E A
A E u E
0 u
(C-2)
because 1
T
u u (assuming that the unit-norm constraint is satisfied throughout the
maneuver to realize physical rotations) and
1 3 1
E u 0 . Also,
190
2 2 2 2
0 3 1 2 0 1 2 3 0 1 2 3
2 2 2 2
0 1 2 3 0 3 1 2 1 3 0 2
4 2 4
2 4 4
2 ,
t
t
t
t
t
t
X X
u u u u u u u u u u u u Y Y
Z Z
X X
u u u u u u u u u u u u Y Y
Z Z
22
2
1
A u
T ξ
T ξ
(C-3)
where in the second equality the relation Eq. (3.16) has been used. Also, it is noted that
when the attitude constraints Eqs. (3.19c) and (3.19d) are satisfied,
22
A u in Eq. (C-3) is
exactly zero.
Likewise,
T
22 1
A E simplifies to
0
2 .
0
3 1 T
22 1
3 2
T ξ T ξ
A E
T ξ T ξ
(C-4)
With Eq. (C-3) and Eq. (C-4), Eq. (C-2) then becomes
2 2 0 2
2 0 2 2 ,
2
2 0 0 0
2 3 1
T
T 22 22 1
22 1 3 2
1 3
T ξ T ξ T ξ
A u A E
A E T ξ T ξ T ξ
0
(C-5)
and the rank of
T
22
A E is three (full rank) if
3
T ξ is not exactly zero. Since Tξ is the
vector connecting the follower satellite and the target expressed along the body frame,
3
T ξ (the z-component of Tξ in the body coordinate frame) is zero only when the z-axis
of the body frame is perpendicular to the vector connecting the follower satellite and the
target (i.e. 90
in Fig. C-1). In typical cases the attitude control begins with a little
error angle between the z-axis of the body and the vector Tξ , and this angle converges
191
to zero by the attitude controller (see Eq. (2.48)), so it is guaranteed that 0
3
T ξ during
the maneuver,
T
22
A E and
22
A has full rank. Since the constraint matrix A has the form of
,
11
21 22
A 0
A
A A
(C-6)
and
22
A has three linearly independent rows, no linear combination of these rows of
22
A
can make any zero rows as required in the (1,2) element of the matrix in the right hand
side of Eq. (C-6). Hence, the matrix A has full rank if
11
A is full rank, and this is the end
of the proof.
Fig. C-1 During attitude control the angle remains small, which guarantees nonzero
3
T ξ .
192
Appendix D
To evaluate the right hand side of Eq. (3.37) the matrix
-1
-1 T
AM A is first computed:
1
1
ˆ
4
1 1
.
1 1 1
ˆ
4
m
m m
m m
T T
3 3 3 4
11 -1 T 11 21
T
T -1 21 22 22
4 3
T T
11 11 11 21
T T T -1 T
21 11 21 21 22 22
I 0
A 0 A A
AM A
A A 0 A
0 E J E
A A A A
A A A A A E J EA
(D-1)
where the matrix E is given in Eq. (2.64), and is orthogonal. The (1,1) component of Eq.
(D-1), which is the s by s matrix
1
m
T
11 11
A A , is invertible if
11
A has full rank, or the
orbital constraints are linearly independent of each other. Let B be the inverse of the
symmetric matrix
-1 T
AM A so that
.
-1
11 12 -1 T
T
12 22
B B
B AM A
B B
(D-2)
Since A has full rank (see Appendix C) and M is positive definite, the matrix
-1 T
AM A
is also positive definite, which guarantees that
1
m
T
11 11
A A (its (1,1) component) and
1 1
ˆ
4 m
T T -1 T
21 21 22 22
A A A E J EA (its (2,2) component) are positive definite, as well. Then,
t
-1
C T -1 T
T T
11 12 1 11 1 11 21
T T
12 22 2 21 22 2 22
T T T T T
11 11 21 12 1 11 1 11 12 21 22 2 21 1 22 2
T T T T
22 12 1 11 1 22 22 2 21 1 22 2
Q A AM A b Aa
B B b A 0 a A A
B B b A A a 0 A
A B A B b A a A B A B b A a A a
A B b A a A B b A a A a
.
(D-3)
193
As noted earlier, the first component (which is a 3-vector) of
t
C
Q is the control force
for orbital control, and the second component (which is a 4-vector) of
t
C
Q is the
generalized quaternion torque for attitude control. To use the result of Eq. (3.37), the
explicit forms of
11
B ,
12
B , and
22
B should be determined first. Hence, the following
well-known result is employed:
The inverse of a partitioned symmetric matrix
,
11 12
T
12 22
V V
V
V V
is the symmetric matrix W (Horn and Johnson [56])
11 12 -1
T
12 22
W W
W V
W W
with
, ,
-1 -1
-1 T -1 -1 T -1 T -1 -1
11 11 12 22 12 11 11 12 22 12 11 12 12 11 12 11 12 22
W V V V V V V V V V V V V V W V V W
.
-1 -1
T -1 -1 -1 T -1 T -1
22 22 12 11 12 22 22 12 11 12 22 12 12 22
W V V V V V V V V V V V V V
By this theorem, the control force 3-vector for orbital control is explicitly determined as
,
t
m m
m
C T T T T T
11 11 21 12 1 11 1 11 12 21 22 2 21 1 22 2
-1 -1
T T T T T
11 11 11 1 11 1 3 3 11 11 11 11 21
-1
-1 -1
T -1 T T T T T T
21 21 22 u 22 21 11 11 11 11 21 2 21 1 22 2 21 11 11 11 1 11 1
F A B A B b A a A B A B b A a A a
A A A b A a I A A A A A
A A A M A A A A A A A b A a A a A A A A b A a
(D-4)
and the generalized quaternion torque 4-vector for attitude control is explicitly
determined as
,
t
m m
C T T T T
u 22 12 1 11 1 22 22 2 21 1 22 2
-1
-1 -1
T T -1 T T T T T T
22 21 21 22 u 22 21 11 11 11 11 21 2 21 1 22 2 21 11 11 11 1 11 1
Γ A B b A a A B b A a A a
A A A A M A A A A A A A b A a A a A A A A b A a
(D-5)
194
where
ˆ
: 4
T
u
M E JE . Eqs. (D-4) and (D-5) are the same as those given in Eqs. (3.38) and
(3.39) upon defining
0
P and
1
P as in Eq. (3.40),
O
α and D as in Eq. (3.41), and
O,A
α as
in Eq. (3.42).
195
Appendix E
Here it is shown that the control force and torque in Eq. (3.46) can be simplified more if
the attitude constraints and the unit-norm constraint, Eqs. (19c), (19d), and (30), are
satisfied. Using Eqs. (3.33)-(3.35), the matrices
21
A ,
22
A , and
2
b are decomposed as
,
21
21
1 3
A
A
0
,
2
22
T 22
A
A
u
,
2N
2
2
b
b
u
(E-1)
where
i
and
i
( 3,4,5 i ) are all zero because the attitude constraints are assumed to
be satisfied. Then, inserting these relations into Eqs. (D-4) and (D-5), one can finally
obtain the following simplified control force and generalized quaternion torque:
,
4
t m m
m
-1
C T T T
11 11 11 1 11 1 1 21
-1
-1
T T -1 T T T
21 1 21 22 1 1 22 2 21 1 1 22 2 21 11 11 11 1
F A A A b A a P A
A P A A E J E A b A P a A a A A A A b
(E-2)
,
4
m
t m
-1
-1
C T T T -1 T T T
u 22 21 1 21 22 1 1 22 2 21 1 1 22 2 21 11 11 11 1
A A P A A E J E A b A P a A a A A A A b
(E-3)
where the fact
22 2 1
A u 0 is used from Eq. (C-3) because the attitude constraints are
satisfied. The physically applied torque t Γ which is a 3-vector is given by Eq. (3.6)
and is obtained as
.
2 4
m m
t
-1
-1
C T T T -1 T T T
1 22 21 1 21 22 1 1 22 2 21 1 1 22 2 21 11 11 11 1
Γ E A A P A A E J E A b A P a A a A A A A b
(E-4)
Next, the following are defined:
, , , t
-1
T T
O 11 11 11 1 11 1
α X X A A A b A a
(E-5)
196
, , , , ,
4
m
t
-1
T T -1 T
t 21 1 21 22 1 1 22
D X X u X A P A A E J E A
(E-6)
and
, , , , , , , . t
O,A t t t 2 21 1 22 2 21 O
α X X u X X u X D b A a A a A α
(E-7)
Using these definitions, Eqs. (E-2) and (E-4) can now be represented in matrix form as
1
.
1
2
t
m
t
T
C
3 3 21
O
C
T
O,A
1 22
I P A
α F
α Γ
0 E A
(E-8)
It must be noted that Eq. (E-8) has been derived on the assumption that the attitude
constraints and the unit-norm constraint are satisfied. If the initial conditions, for
example, do not meet the attitude constraints (i.e. the z-axis of the body frame does not
point to a desired target at the initial time), a more general control force and torque (Eq.
(3.46)) should be used.
Abstract (if available)
Abstract
This dissertation handles two topics that are fundamental in understanding aerospace and mechanical systems. In the first topic, a new controller for maintaining a desired formation of multiple satellites is developed in the presence of attitude tracking requirements under various perturbations and model uncertainties. Both orbital and attitude dynamics are simultaneously considered and continuous thrust propulsion systems are assumed. Unlike most studies, no linearizations and/or approximations are made in dynamics or the controllers. The novel controller developed herein provides a remarkable improvement over the current state-of-theart in that the control force and torque to satisfy the given orbital and/or attitude constraints are obtained in completely closed form. This new, analytical solution can be easily used for on-orbit, real-time control with low computational burden. Furthermore, it is useful in estimating the magnitude of the required control inputs and results in some interesting consequences, including the separation principle, which describes when the orbital control can be completely decoupled from the attitudinal control required to be applied to the follower satellite. Also, an additional additive controller that compensates for model uncertainties is developed. This is done by using the desired trajectory of the nominal system as the tracking signal, and is based on a generalization of the concept of sliding surfaces. The resulting closed-form control causes the desired attitude and orbital requirements of the nominal system to be met in the presence of unknown, but bounded, model uncertainties. ❧ The second topic deals with the inverse problem of the calculus variations. To develop the equations of motion for a mechanical system, it is usual to define a Lagrangian function first and then obtain the equations of motion by applying the Lagrange equation when the forces are derivable from a potential function. However, it is interesting and useful to investigate the types of forces that can engendered when substituting the forces into the Lagrange equation guarantees the proper equations of motion even when the forces do not arise from a potential, that is, when they are nonconservative. By extending the existing Bolza’s approach, a more general and systematic way of handling the problem is found for single-degree-of-freedom systems. Then, new results for linearly damped multi-degree-of-freedom systems are obtained using extensions of results for single-degree-of-freedom systems. Conservation laws for these damped multi-degree-of-freedom systems are also found using the Lagrangians obtained. These new findings can be used for development of approximate solutions of nonlinear differential equations or various numerical techniques.
Linked assets
University of Southern California Dissertations and Theses
Conceptually similar
PDF
Flexible formation configuration for terrain following flight: formation keeping constraints
PDF
New approaches in modeling and control of dynamical systems
PDF
Optimal guidance trajectories for proximity maneuvering and close approach with a tumbling resident space object under high fidelity J₂ and quadratic drag perturbation model
PDF
An analytical dynamics approach to the control of mechanical systems
PDF
Techniques for analysis and design of temporary capture and resonant motion in astrodynamics
PDF
Transient modeling, dynamic analysis, and feedback control of the Inductrack Maglev system
PDF
Application of the fundamental equation to celestial mechanics and astrodynamics
PDF
On the synthesis of controls for general nonlinear constrained mechanical systems
PDF
On the synthesis of dynamics and control for complex multibody systems
PDF
Control of spacecraft with flexible structures using pulse-modulated thrusters
PDF
Dynamic modeling and simulations of rigid-flexible coupled systems using quaternion dynamics
PDF
Trajectory mission design and navigation for a space weather forecast
PDF
Mathematical model and numerical analysis of a shape memory polymer composite cantilever beam for space applications
PDF
Controlled and uncontrolled motion in the circular, restricted three-body problem: dynamically-natural spacecraft formations
PDF
Terrain-following motion of an autonomous agent as means of motion planning in the unknown environment
PDF
The development of an autonomous subsystem reconfiguration algorithm for the guidance, navigation, and control of aggregated multi-satellite systems
PDF
Computationally efficient design of optimal strategies for passive and semiactive damping devices in smart structures
PDF
Development and evaluation of high‐order methods in computational fluid dynamics
PDF
Optimal design, nonlinear analysis and shape control of deployable mesh reflectors
PDF
Solution of inverse scattering problems via hybrid global and local optimization
Asset Metadata
Creator
Cho, Hancheol
(author)
Core Title
New approaches to satellite formation-keeping and the inverse problem of the calculus of variations
School
Viterbi School of Engineering
Degree
Doctor of Philosophy
Degree Program
Aerospace Engineering
Publication Date
10/19/2012
Defense Date
07/06/2012
Publisher
University of Southern California
(original),
University of Southern California. Libraries
(digital)
Tag
formation flying,inverse problem of the calculus of variations,Lagrangian,OAI-PMH Harvest,orbital/attitude control,satellite
Language
English
Contributor
Electronically uploaded by the author
(provenance)
Advisor
Udwadia, Firdaus E. (
committee chair
), Erwin, Daniel A. (
committee member
), Yang, Bingen (Ben) (
committee member
)
Creator Email
hancheoc@usc.edu,hancheol.cho@gmail.com
Permanent Link (DOI)
https://doi.org/10.25549/usctheses-c3-106916
Unique identifier
UC11290220
Identifier
usctheses-c3-106916 (legacy record id)
Legacy Identifier
etd-ChoHancheo-1257.pdf
Dmrecord
106916
Document Type
Dissertation
Rights
Cho, Hancheol
Type
texts
Source
University of Southern California
(contributing entity),
University of Southern California Dissertations and Theses
(collection)
Access Conditions
The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...
Repository Name
University of Southern California Digital Library
Repository Location
USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA
Tags
formation flying
inverse problem of the calculus of variations
Lagrangian
orbital/attitude control
satellite