<%BANNER%>

Dynamics and Control for Formation Flying System


PAGE 1

DYNAMICS AND CONTROL FOR FORMATION FLYING SYSTEM By YUNJUN XU A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2003

PAGE 2

Copyright 2003 by Yunjun Xu

PAGE 3

Dedicated to my parents, who always loved; my wife, Yanli who always encouraged; my son, Ronald who always delighted.

PAGE 4

ACKNOWLEDGMENTS I would like to express my gratitude to my advisor, Dr. Norman Fitz-Coy, for his support, patience, and encouragement throughout my graduate studies. It is not often that one finds an advisor and colleague that always find the time to listen to the little problems and roadblocks that unavoidably crop up in the course of performing research. His technical and editorial advice was essential to the completion of this dissertation and has taught me innumerable lessons and insights on the workings of academic research in general. My thanks also go to the members of my committee, Drs. Andrew Kurdila, Rick Lind, Peter Ifju, and Michael Nechyba, for reading previous drafts of this dissertation and providing many valuable comments that improved the presentation and contents of this dissertation. I would like to thank Dr. Paul Mason for his help during my whole research period. The friendship of colleagues in the Dynamics and Control Group is much appreciated and has led to many interesting and good-spirited discussions relating to this research. Last, but not least, I would like to thank my wife, Yanli, and my son, Ronald, for their understanding and love during the past few years. Their support and encouragement were in the end what made this dissertation possible. My parents receive my deepest gratitude and love for their dedication and the many years of support that provided the foundation for this work. iv

PAGE 5

TABLE OF CONTENTS page ACKNOWLEDGMENTS..iv LIST OF TABLESviii LIST OF FIGURES...ix ABSTRACT...xvii CHAPTER 1 INTRODUCTION....1 1.1 Motivation.1 1.2 Research Scope.8 1.3 Dissertation Outline..9 2 MULTIPLE VEHICLES DYNAMIC MODEL..11 2.1 Introduction..11 2.2 Discussion Of Multiple Vehicles Dynamic Model 2.3 General 2-Body Relative Motion Models..17 2.3.1 Relative Motion Dynamic Model (Point Mass Case)... ...18 2.3.2 Simulation Results and Discussion (Point Mass Case)..28 2.3.2.1 Initial condition.....28 2.3.2.2 Results and discussion.....28 2.3.3 Specialized Relative Motion Model (Point Mass Case) 2.3.3.1 Linear model (Clohessy-Wiltshire equation).34 2.3.3.2 Circular reference orbit nonlinear model...36 2.3.3.3 Elliptic reference orbit nonlinear model.37 2.3.4 Simulation Results and Discussion (Point Mass Case)...38 2.3.5 Relative Motion Dynamic Model (Rigid Body Case)....40 2.3.5.1 Desired attitude..40 2.3.5.2 Attitude model....47 2.3.6 Simulation Results and Discussion (Rigid Body Case)..49 2.4 Perturbation 2.4.1 Atmospheric Drag..52 2.4.2 Earth Oblateness.53 v

PAGE 6

2.4.3 Perturbation Effects on the Generalized Dynamics Model 2.4.4 Perturbation Effects on Keplerian Elements..58 2.5 Conclusion...64 3 CONTROL STRATEGIES: ORBIT..65 3.1 Introduction 3.2 Performance Criterion Of The Leader-Follower Pair 3.3 Relative Position Controller Design...70 3.3.1 Linear-Quadratic State Feedback Regulator (LQR)..70 3.3.1.1 LQR 3.3.1.2 LQR synthesis model.71 3.3.1.3 LQR analysis model...74 3.3.2 /H Controller..........76 3.3.2.1 Structured singular value ( ) theory 6 3.3.2.2 synthesis model..77 3.3.2.3 analysis model81 3.3.2.4 NP, RS, and RP..81 3.3.3 Feedback Linearization..83 3.3.3.1 Feedback linearization....83 3.3.3.2 Feedback linearization analysis model89 3.3.4 Sliding Mode Controller (SMC)0 3.4 Relative Position Simulation Results and Discussion....92 3.4.1 State Feedback Gain 93 3.4.2 Robustness and Performance Analysis..95 3.4.3 Time Response Based on C-W Model96 3.4.4 Time Response Based on C-W Model with Uncertainties.....98 3.4.5 Time Response Based on Actual Nonlinear Model with Thruster Limitation.101 3.4.6 Conclusion2 4 CONTROL STRATEGIES: ATTITUDE AND ORBIT.104 4.1 Introduction..104 4.2 Relative Attitude Control.104 4.2.1 Control Architecture 4.2.2 The Inner Loop Controller106 4.3 Combined Relative Position and Attitude Controller .112 4.4 Boundary Layer Control..116 vi

PAGE 7

vii 4.4.1 Hard Switching Control117 4.4.2 Boundary Layer Control..117 4.4.3 Simulation Results...118 4.5 Genetic Algorithm125 4.5.1 Optimization Problem..125 4.5.2 Genetic Algorithm127 4.5.3 Simulation Results131 5 RECONFIGURATION CONTROL LER IN FORMATION FLYING..136 5.1 Introduction.136 5.2 Formation Reconfiguration Strategy 5.3 Perceptive Frame Control (PFC).....143 5.4 Controller Design.145 5.4.1 Perceptive Frame Controller145 5.4.2 Formation Maintenance Controller.146 5.4.3 Formation Reconfiguration Controller147 5.4.4 SIMULINK Block Diagram 5.5 Results .149 5.6 Summary.156 6 CONCLUSION AND FUTURE WORK 6.1 Conclusion157 6.2 Future Work APPENDIX QUATERNION AND LIE DERIVATIVE.160 A.1 Quaternion.......160 A.2 Lie Derivative..161 LIST OF REFERENCES.162 BIOGRAPHICAL SKETCH170

PAGE 8

LIST OF TABLES Table page 2-1 Different cases been simulated 2-2 Classical elements for an Earth satellite..59 3-1 Gain and phase margins...95 3-2 Robustness performance and controller order comparison..96 3-3 Fuel consumption.98 3-4 Fuel consumption...101 4-1 Cases..107 4-2 Cases..110 4-3 Cases..113 4-4 GA parameters definition and the data used 4-5 Simulation results of GA + SMC + BLC..133 5-1 The initial orbit information of the formation satellites 5-2 Orbit information of R1 and R2 in figures 3 and 4 5-3 Formation maintenance controllers for follower 1 and follower 2 5-4 Formation reconfiguration controllers for redundant 1 and redundant 2..147 5-5 Configuration code for in-plane formation flight simulation viii

PAGE 9

LIST OF FIGURES Figure page 1.1 Small satellite for DS1....................................................................................................4 1.2 The probe and mission outline........................................................................................5 1.3 The satellite for EO-1......................................................................................................6 1.4 Formation flying of ST5.................................................................................................7 1.5 Size of the micro-satellite and the communication architecture.....................................7 2.1 Classification of the relative dynamics model................................................................11 2.2 N bodies system..............................................................................................................15 2.3 Zones of variations..........................................................................................................16 2.4 Relative dynamic model.................................................................................................18 2.5 Coordinate system of the relative position model...........................................................19 2.6 Coordinate system of the leader satellites frame...........................................................19 2.7 Orbital elements..............................................................................................................25 2.8 Orbital information.........................................................................................................31 2.9 Relative position in x direction.......................................................................................31 2.10 Relative position in y direction.....................................................................................31 2.11 Relative position in z direction.....................................................................................31 2.12 Relative position in 3D.................................................................................................32 2.13 Model difference...........................................................................................................32 2.14 Case 9 model difference...............................................................................................32 2.15 Case 20 model difference..............................................................................................34 ix

PAGE 10

2.16 Relative position in x direction.....................................................................................39 2.17 Relative position in y direction.....................................................................................39 2.18 Difference between CW and generalized model..........................................................39 2.19 Satellite attitude body frame.........................................................................................41 2.20 Spherical right triangle..................................................................................................42 2.21 Desired attitude I...........................................................................................................43 2.22 Desired attitude II.........................................................................................................43 2.23 Desired attitude III........................................................................................................44 2.24 Error angles in the coordinate system. EL is the elevator error angle, whereas, AZ is the azimuth error angle.....................................................................................46 2.25 Theta 1 (EL) and Theta2 (AZ) in radians.....................................................................50 2.26 Theta 1 (EL) and Theta2 (AZ) in radians.....................................................................50 2.27 Theta1 (EL) and Theta2 (AZ).......................................................................................51 2.28 Atmosphere drag...........................................................................................................53 2.29 The perturbation stems from the Earth oblateness........................................................54 2.30 Atmosphere drag effects...............................................................................................57 2.31 J 2 perturbation effects...................................................................................................57 2.32 Combined perturbation effects......................................................................................57 2.33 Error in position............................................................................................................58 2.34 Error in attitude.............................................................................................................58 3.1 Combined state space......................................................................................................72 3.2 Analysis model of LQR controller without uncertainty.................................................75 3.3 Analysis model of LQR controller with uncertainty.......................................................75 x

PAGE 11

3.4 Typical loop gain and phase margin model....................................................................76 3.5 Open loop model with parametric uncertainty................................................................79 3.6 Synthesis Model..............................................................................................................80 3.7 Synthesis model input-output relation............................................................................81 3.8 Analysis model for nominal system................................................................................81 3.9 Analysis model for actual system...................................................................................81 3.10 Closed-loop model with uncertainty.............................................................................82 3.11 Equivalent model 1.......................................................................................................82 3.12 Equivalent model 2.......................................................................................................82 3.13 Closed-loop analysis model based on the feedback linearization.................................89 3.14 Closed-loop analysis model for the linear system........................................................90 3.15 Relative position error in x...........................................................................................96 3.16 Relative position error in x...........................................................................................96 3.17 Relative position error in y...........................................................................................97 3.18 Relative position error in y...........................................................................................97 3.19 Relative position error in z............................................................................................97 3.20 Relative position error in z............................................................................................97 3.21 Control force in x..........................................................................................................98 3.22 Control force in y..........................................................................................................98 3.23 Control force in z..........................................................................................................98 3.24 Relative position error in x...........................................................................................99 3.25 Relative position error in x...........................................................................................99 xi

PAGE 12

3.26 Relative position error in y...........................................................................................99 3.27 Relative position error in y...........................................................................................99 3.28 Relative position error in z............................................................................................100 3.29 Relative position error in z............................................................................................100 3.30 Control force in x..........................................................................................................100 3.31 Control force in y..........................................................................................................100 3.32 Control force in z..........................................................................................................101 3.33 Relative position error in x...........................................................................................101 3.34 Relative position error in y...........................................................................................101 3.35 Relative position error in z............................................................................................102 4.1 Attitude control architecture...........................................................................................105 4.2 Relative position difference between exact and uncertainties model.............................108 4.3 Relative position difference between exact and uncertainties model.............................108 4.4 Elevator error angle difference between the exact and uncertainties models.................109 4.5 Azimuth error angle difference between the exact and uncertainties models................109 4.6 Elevator error angle if no control is applied...................................................................109 4.7 Azimuth error angle if no control is applied...................................................................109 4.8 Elevator error angle under the sliding mode controller..................................................111 4.9 Azimuth error angle under the sliding mode controller in transient stage from 0 ~ 600 seconds.................................................................................................................111 4.10 Azimuth error angle under the sliding mode controller from 5000 ~ 5800 seconds....111 4.11 Elevator error angle in case 5........................................................................................111 4.12 Azimuth error angle in case 5.......................................................................................112 xii

PAGE 13

4.13 Relative position in x direction.....................................................................................113 4.14 Relative position in y direction.....................................................................................113 4.15 Relative position in z direction.....................................................................................114 4.16 The elevation error angle..............................................................................................114 4.17 Zoomed in elevation error angle...................................................................................114 4.18 Zoomed in azimuth error angle.....................................................................................114 4.19 Zoomed in azimuth error angle.....................................................................................114 4.20 Force commands in x direction.....................................................................................114 4.21 Zoomed in force commands in x direction...................................................................115 4.22 Force commands in y direction.....................................................................................115 4.23 Zoomed in force commands in z direction...................................................................115 4.24 Torque commands in axis 1..........................................................................................115 4.25 Torque commands in axis 2..........................................................................................116 4.26 Sliding surface in SMC ................................................................................................118 4.27 Sliding surface in SMC + BLC.....................................................................................118 4.28 Relative position in x direction, transient stage............................................................119 4.29 Relative position error in x direction, stable stage........................................................119 4.30 Relative position in y direction, transient stage............................................................119 4.31 Relative position error in y direction, stable stage........................................................119 4.32 Relative position in z direction, transient stage............................................................119 4.33 Relative position in z direction, stable stage.................................................................119 4.34 Sliding surface for the relative position in x direction, transient stage.........................120 xiii

PAGE 14

4.35 Sliding surface for the relative position in x direction, stable stage.............................120 4.36 Sliding surface for the relative position in y direction, transient stage.........................120 4.37 Sliding surface for the relative position in y direction, stable stage.............................120 4.38 Sliding surface for the relative position in z direction, transient stage.........................120 4.39 Sliding surface for the relative position in z direction, stable stage.............................120 4.40 The elevator error angle in the transient stage..............................................................122 4.41 The elevator error angle in the stable stage..................................................................122 4.42 The azimuth error angle in the transient stage..............................................................122 4.43 The azimuth error angle in the stable stage..................................................................122 4.44 Force commands in x direction of the local frame, in transient stage..........................123 4.45 Force commands in x direction of the local frame, in stable stage...............................123 4.46 Force commands in y direction of the local frame, in transient stage..........................123 4.47 Force commands in y direction of the local frame, in stable stage...............................123 4.48 Force commands in z direction of the local frame, in transient stage...........................124 4.49 Force commands in z direction of the local frame, in stable stage...............................124 4.50 Torque commands in body axis 1, in transient stage....................................................124 4.51 Torque commands in body axis 1, in stable stage........................................................124 4.52 Torque commands in body axis 2, in transient stage....................................................124 4.53 Torque commands in body axis 2, in stable stage........................................................124 4.54 A generation of the chromosome..................................................................................128 4.55 Evaluation, fitness, and selection..................................................................................129 4.56 Crossover......................................................................................................................130 xiv

PAGE 15

4.57 Mutation........................................................................................................................130 4.58 The relative position error in x-direction......................................................................134 4.59 Control force in x-direction..........................................................................................134 4.60 The relative position error in y-direction......................................................................134 4.61 Control force in y-direction..........................................................................................134 4.62 The relative position error in z-direction......................................................................134 4.63 Control force in z-direction...........................................................................................134 4.64 Elevation error angle.....................................................................................................135 4.65 Torque in body axis 1...................................................................................................135 4.66 Azimuth error angle......................................................................................................135 4.67 Torque in body axis 2...................................................................................................135 5.1 Absolute..........................................................................................................................137 5.2 Leader-follower Structure...............................................................................................137 5.3 Absolute and interconnection.........................................................................................137 5.4 Sub formation..................................................................................................................137 5.5 3-satellite formation architecture with two redundant satellites.....................................139 5.6 If F 2 fails.........................................................................................................................141 5.7 If F 1 fails.........................................................................................................................142 5.8 Perceptive frame controllers structure...........................................................................145 5.9 SIMULINK program diagram........................................................................................148 5.10 Relative position in x local frame. (a) whole range, and (b) from 5670 7500 second..................................................................................................................150 xv

PAGE 16

xvi 5.11 Relative position in y local frame, (a) whole range, and (b) from 5500-10000 seconds.................................................................................................................150 5.12 Relative position differences in z lo cal frame, (a) 0-800 seconds, and (b) 5650-6500 seconds.................................................................................................................151 5.13 Force commands in x local frame. (a) whole range, and (b) 5600-10000 seconds......152 5.14 Force commands in y local frame. (a) whole range, and (b) 5600-9000 seconds........152 5.15 Force commands in z local frame. (a) 0-800 seconds, and (b) 5900-6300 seconds.....153 5.16 Elevation error angle. (a) w hole range, and (b) 5670-15000 seconds..........................154 5.17 Azimuth error angle. (a) w hole range, and (b) 5900-7000 seconds..............................154

PAGE 17

Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy DYNAMICS AND CONTROL FOR FORMATION FLYING SYSTEM By Yunjun Xu December 2003 Chair: Norman Fitz-Coy Major Department: Mechanical and Aerospace Engineering Due to the potential advantages of the micro-satellites formation flying system, such as cost reduction, flexibility, observational baseline increase, and reliability improvement, numerous formation related missions have been proposed and seriously studied. Along these many potential benefits, come a variety of challenges. Dynamics and control of the multi-satellites in the formation is one of these issues and is addressed in this dissertation. To date, most formation flying dynamics and control studies are based on the linear Clohessey-Wiltshire equations and none address the problem with a coupling between the attitude and orbit motions. Furthermore, commercial software, such as STK and FreeFlyer, is not capable of feedback control simulation. These observations motivate this research to develop a generalized coupled attitude-orbit relative motion and to design controllers based on this model. The developed generalized dynamics model is based on the leader-follower structure and the desired attitude is coupled with the real time relative position due to the formation mission requirement. This nonlinear model is applicable to any closed xvii

PAGE 18

reference orbit cases and incorporates a coupling between the attitude and orbit. Also, atmospheric drag and J 2 perturbations are included. The formation flying control system can be regarded as a two-tiered architecture. In the lower level, continuous controllers are designed to meet the objectives of each vehicle, and in the higher level, discrete time event determines the formation configuration and reconfiguration. After comparison with LQR, /H synthesis, and feedback linearization strategies, sliding mode control is chosen as the basic control method for the lower level in this dissertation. Furthermore, boundary layer control and genetic algorithm are implemented in order to reduce chatters with the implementation of the sliding mode control and reduce the fuel consumption. In the upper level, perceptive frame concept is adopted for the discrete event determination. It is demonstrated that the proposed control architecture is capable of tracking and maintaining the desired position and attitude under different formation missions. xviii

PAGE 19

CHAPTER 1 INTRODUCTION 1.1 Motivation In the 1950s, the first satellite was launched by the Soviet Union [1]. From that moment on, more and more satellites have been launched into space for the purpose of Earth observation, deep space exploring, military surveillance, navigation (GPS), commercial and military communication, and weather prediction. To date, conventional monolithic satellites typically have the following deficiencies. The first deficiency is their large size and/or weight. A significant portion of a mission's budget is spent on a launch vehicle. The larger a satellite, the larger and more costly a launch vehicle is required to boost it into space. The second deficiency of the conventional satellite is its lack of flexibility for missions. When the satellite is located in its orbit, the mission is fixed. The instruments of the satellite will limit the observing baseline. The last deficiency is the bad redundancy. Whenever a failure occurs on the satellite, typically the entire mission fails or the redundant subsystem may need to be used, which will increase the payload of the satellite. Therefore, a new concept, formation flying system or distributed space system, was introduced in the early 80s by Hartman and Hadaegh et al. [2, 3]. In 1998, the formation flying system has received renewed interest by numerous government agencies. The formation flying systems are multiple mini, micro, or nano satellites that can work together in order to perform a unified mission [4]. Such a system will take the 1

PAGE 20

2 place of the conventional mono lithic satellites. This system can range from global constellations offering extended service coverage to clusters of highly coordinated vehicles that perform distributed sensing. [4]. NASAs New Millennium Program [5] propos es the properties of the formation flying system as the following: Fleets of vehicles Vehicles interact and cooperate with one another to achieve mission goals Vehicle pointing managed collectively Vehicle positioning managed collectively Fleets evolve over time to extend and expand capability Self-controlling vehicles eliminate the need for extensive ground support With respect to the conventional single large satellites, a formation flying system may provide numerous advantages [5-8]. So me of these advantages are summarized below. Cost Reduction: Due to the size and wei ght reduction, the launch cost will decrease substantially. Using the Space Shuttle as an ex ample, the cost to deploy a payload into a low-Earth orbit is about $13,200 per kilogr am ($6,000 per pound) [5]. Most current satellites have masses on the order of 1,000 K g. But in formation flying system, the mass of a micro-sat will be about 100 Kg, and that of a nano-sat will be approximately 20 Kg. Therefore, at least 50 nano-satellites can be la unched into the space at the same cost with potentially more powerful functions. More flexible, efficient use of resources: The formation of the multiple satellites can be changed autonomously or manually according to different missions or error conditions. Extensive co-observing programs can be conducted without using extensive ground support. For instance, in the leader-follower communication architectu re, the satellites do not need to communicate with th e ground station all the time. On ly the leader satellite is in charge of the task and the autonomous handling software will determine when the communication with the ground support is necessary.

PAGE 21

3 Increase precision and observa tional baselines. For exampl e, A team of powerful XRay telescopes unlike anything this part of the universe ha s ever seen. These telescopes will work in unison to simultaneously observe the same distant objects, combining their data and becoming 100 times more powerful th an any single X-ray telescope that has come before it.[5] Enhanced survivability and increased reliability: In a case in which a certain number of the satellites in the formation system faile d, the mission may still be accomplished. This is an important advantage over the conventional monolithic satellite. With these potential benefits, however, co me a variety of t echnical challenges. Some of them are the following: 1. Almost all space missions are flown with a single spacecraft. Controlling spacecraft or satellites in formation flight is very complicated. 2. Hybrid system and reconfi guration control is required. 3. The ground-to-multiple-sate llite communications path is more complex than the ground-to-single-satellite communications. 4. Because of the limitation in size and we ight, designing and manufacturing the mini, micro, or nano size components for the satellit es in formation flying system are not easy issues [9]. Although complete formation flying systems have not been implemented in space, NASA and other government agen cies have begun to launch several missions in order to validate critical technical issues, and to incr ease the confidence of this brand new space exploration era. NASAs New Millennium Program, which includes such missions as Deep Space 1, Deep Space 2, Earth Obse rving 1, Earth Observing 3, and Space Technology 5 [5, 10-12], is an example of a government program that is specific to the spacecraft formation. These mi ssions will be briefly discussed. The Deep Space 1 (DS1) mission was launched on October 24, 1998. The retired date of DS1 was December 18, 2001. Figure 1.1 shows the picture of the small satellite. The total mass of DS1 is 486.3 Kg. Several te chniques for the multiple satellites have

PAGE 22

4 been tested in this mission; for example, low-mass communications devices, which are half as heavy and 1/3 as costly as before, microelectronics and spacecraft structure, which tried to demonstrate futuristic technologies for making the spacecraft smaller, lighter, and more efficient, and autonomous operations system that can plan, make decisions, and operate by itself. Figure 1. 1 Small satellite for DS1 The Deep Space 2 (DS2) project includes two probes, piggybacked aboard the Mars Polar Lander. It was launched on January 3, 1998. The two probes were slammed into Martian soil on December 3, 1999. Figure 1.2 shows the probe and the mission diagram. DS2 is aimed to design, develop, and flight validate advanced systems and instruments in deep space in order to lower the cost and risk for future formation missions. The total mass of the probe is 6.5 Kg [13].

PAGE 23

5 Figure 1. 2 The probe and mission outline Although it is in still in an elementary stage, EO-1 is actually the first prototype in the New Millennium Program of a formation flying system. EO-1 works with the LandSat-7 to form a simple formation once or twice each day. Landsat 7 is a U.S. satellite used to acquire remotely sensed images of the Earths land surface and coastal regions. At this point, EO-1 and LandSat-7 will image the same point in order to get the piecewise continuous data. The EO-1 small satellites orbit is a 705 Km circular, sun-synchronous orbit inclined at a 98.7 0 to the equator. It is still in an elementary stage of the formation-flying concept because EO-1 can only coordinate satellites when they are in one or two points of the orbit and the complicated technologies, such as multiple vehicles communication, fault tolerant algorithm, and autonomous system, have not been tested. The mission started on November 21, 2000. Figure 1.3 shows a picture of the EO-1s small satellite. The mass of the small satellite is 90 Kg.

PAGE 24

6 Figure 1. 3 The satellite for EO-1 The objective of the EO-3, which is also called GIFT-IOMI [14], is to improve the observation of the temperature, moisture, wind velocity, and test the new remote sensor, which will be used in future missions. The planned launch date is 2005~2006. The mass of the instrument is about 100 Kg. Satellite Technology 5 (ST-5) is composed of three miniature satellites with a mass of 19.4 Kg and 54.2 cm X 28.6 cm dimensions. Although it is small, it has full service including guidance, navigation, control, propulsion, bandwidth, and communication systems. This mission will be launched in 2004.

PAGE 25

7 Figure 1. 4 Formation flying of ST5 Figure 1. 5 Size of the micro-satellite and the communication architecture. The right-hand side of Figure 1.5 shows the communication architecture that will be used in this mission. Every individual satellite will communicate with the ground base. In addition to the NASA related projects, there are others working in this area, for example, Constellation X-ray, ACE [15], UoSAT-12, LISA and MAXIM. Furthermore, there are some other projects related to the formation flight systems in other application, such as TechSat 21 [16] of DoD and Air Force formation flight system [17].

PAGE 26

8 A formation flying system is a complicated comprehensive system with numerous sub-systems. The ORION micro-satellite is a correct example of an element of a formation flying system. Its power, communication, dynamics and control chip, GPS receiver, structure, orbit r oute, orbit condition, shape of the formation, and sensor are discussed in Ref. [18]. In this dissertation, it is impossible to cover all these topics. In the next section, the specific research topics of a formation flying system addressed in this dissertation will be expressed. 1.2 Research Scope In this section, the specific research topics that are cove red in this dissertation will be summarized, and a brief introduction of each is presented. A detailed literature review of each topic is presented in the individual chapter. Although the component technologies may be ready, the management and control of a distributed constellati on is still a major area of co ncern. Therefore, the dynamic models and the controllers are two of the key issues discussed in th is dissertation. The literature abounds with papers that use the Clohessey-Wiltshire (C-W) equations as the dynamic model upon which the controller is de signed [19-30]. C-W equation, however, is linear relative motion model that is restricted to the near circular orbi ts, and thus neglects a large class of orbit cases. Recently, some papers [30, 31] studied the nonlinear relative motion problem, but are still restricted to the circular reference orbit case. In this dissertation, the dynamic models are extende d to the generalized, non-circular, closed orbits and are not linearized. These models include gravitational perturbations due to oblateness of the central body as well as atmospheric drag.

PAGE 27

9 The design of the controllers for the le ader-follower hierarchical structure of satellite formation is also a key technology which needs to be solved. To date, most designs are based on the PID or LQR cont roller methodology. Recently, some elements of robust control have been utilized in the controller design [4, 32-34]. In this thesis, nonlinear controllers based on sliding mode c ontrol theory will be utilized. Various controller designs will be evaluated based on performance criteria that address the fuel consumption (proportional to magnitude of co ntrol effort) and the errors between the actual and desired states. Due to planed mission or unplanned changes (e.g., satellite failu re), the formation may need to be reconfigured. A hybrid control architecture will be used to address the reconfiguration issue. Hybrid control me thodology includes a multi-tiered control architecture that has been us ed in the smart highway system [35]. In Stadter [36], the high-level tier uses a discre te even command to achieve formation reconfiguration. Graphs, matrix inequalities, a nd switching are also methods to deal with the system level problem [37]. In this dissertation, a perceptiv e frame concept will be utilized and it will be shown to have a good and smooth performance [38, 39]. 1.3 Dissertation Outline This dissertation is organized as follows. In Chapter 1, the basic concept and the advantages of the satellite formation flying sy stem were addressed. The state-of-the-art of the formation-flying projects was illustrated. Also the research scope of this dissertation will be provided. Chapter 2 addresses the detailed dynamic models needed to analyze a spacecraft formation flying system. The existing models in the literature will be presented and discussed; their deficiencies will be high lighted. After that, the generalized nonlinear

PAGE 28

10 dynamic model of the leader-follower satellites pair will be derived from the first principles. It will be shown that the specialized models that appeared in the literature can be derived from the generalized model under proper conditions. Also, in this chapter, the combined attitude and orbit model will be developed. Finally in Chapter 2, the perturbation terms and their effects on the classical six-orbital elements will be discussed. Chapters 3, 4, and 5 address the control strategies for the formation flying system. A two-tiered control structure is utilized in the formation system. The lower tier addresses the specific vehicle control actions and the upper tier addresses the hierarchy and reconfiguration control actions. The Chapters 3 and 4 will be focused on designing the lower part, which is the controller for the relative position and attitudes maintaining and tracking for the leader and follower satellites pair. The controller study will be based on sliding mode control methodology, and the performance will be compared with the LQR, synthesis, and feedback linarization methods, which often appeared in the literature. All controllers designed are rigid spacecraft approach. Furthermore, some improvements for the sliding mode control are discussed. Chapter 5 addresses the higher tier of the controller and focused on reconfiguration of the whole formation flying system. The controller design is based on the perceptive frame concept. It will be shown that the overall system is stabilized and controllable. Chapter 6 provides some conclusions of this research as well as the suggestions for future work.

PAGE 29

CHAPTER 2 MULTIPLE VEHICLES DYNAMIC MODEL 2.1 Introduction In order to effectively analyze orbit maintenance and re-configuration control strategies for the satellite formation systems, mathematical models of the relative motions between satellites are required. In the past two decades, numerous models have appeared in the literature, with essentially the same assumptions. Typically, a leader/follower architecture is chosen and the leaders orbit information is known [2, 21, 40, 41]. These models can be roughly classified as shown in Figure 2.1. Figure 2. 1 Classification of the relative dynamics model The most widely used model for the relative dynamics motion in the formation flying system is the classical Clohessey-Wiltshire (C-W) equation [19-24, 26-30, 42]. Vassar and Sherwood [21], Redding and Adams [22], and Sedwich et al. [23] were the 11

PAGE 30

12 first to apply this linear model to formation flight systems. Taking advantage of the uncoupling between the in-plane and out-of-plane motions, a closed form controller was designed. In 1989, Yan et al. [26] made some modifications to this linear model. One is discretizing the continuous model, and the other is using v (impulse) instead of normal forces. Through these changes, the C-W solution is represented as 1ii x xv where and i represent discrete time steps, i 1 is the state transition matrix, and is the input matrix. Nevertheless, the results obtained from this modified C-W model are not acceptable for precision formation flight. For example, Vassar and Sherwood [21] stated that the precision of the relative position was on the order of Km. This precision is too large for the types of formation flight missions being proposed. Furthermore, this model is limited to the cases where the leader satellite is in a circular orbit and the disturbance varies slowly with respect to time. 21 Queiroz et al. and Kapila et al. [31, 32, 43] used a nonlinear system as the relative position model, which is also restricted to the condition that the leader satellites orbit is a circular orbit. Also imposed are the conditions that the mass of the satellite and the disturbance experienced by the satellite vary slowly with respect to time. Lyapunov-typed controllers are designed for these nonlinear models and the control precision is on the order of 1 m. Inalhan, Tillerson, and How [44, 45] in January 2002 described the relative position motion for the case of an elliptic reference orbit, in place of the circular reference orbit. In the noncircular case, the time dependence was replaced with the true anomaly dependence. Inalhans paper focuses on developing T-periodic solutions, which refer to having the vehicles return to the initial relative states at the end of one orbital

PAGE 31

13 period. Thus, the dynamic model is not general enough to study the formation system. First it is not applicable to the singular orbit configurations. Singularity occurs when either the inclination or the eccentricity of the leader or followers orbit becomes zero. Second the method only guarantees the two boundary points to satisfy the given condition. Another approach to evaluate elliptic orbit reference case was provided by Schiff et al. [46]. This approach, however, is limited to nearly circular orbits since it is based on the C-W equation, which is implemented in AI Solution, Inc.s software FreeFlyer. This model fails for orbits with any appreciable eccentricity. All of the studies mentioned above consider the vehicles as point masses, and as such, only consider the formation from a relative position point of view. Ahmed et al. [47], Wu and Chen [48], Kim et al. [49], and Xing and Parvez [50] studied the rigid body cases and represented the attitude model using JJ M The desired state of the attitude system is expressed as a function of the angular velocity. The basic idea is to maintain or track a desired angular velocity profile. Based on this kind of model, different control methods were proposed including adaptive asymptotic tracking, /, and sliding mode control. To date, there have been no publications demonstrating a coupled attitude/orbit problem. 2H H Satellite Tool Kit (STK) from Analytical Graphics, Inc., provides the commercial analytical tool for more than 30,000 aerospace, defense, and intelligence professionals. It gives very good vision and open-loop simulation environments for spacecraft-related research and design. However, the closed loop or control input to the spacecraft cannot be implemented up to the most recent version. Therefore, the coupled dynamics model is developed and addressed in this dissertation.

PAGE 32

14 In Section 2.2 of this chapter, the N-body problem is simplified to the 2-body dynamics model. Then the generalized relative dynamics model, both point mass case and rigid body case, are derived in Section 2.3. The special cases, including all the existing models in the literature, are also described in Section 2.3. Finally, perturbation and their effects on the classical orbital elements are discussed in Section 2.4. 2.2 Discussion Of Multiple Vehicles Dynamic Model Satellites experiencing formation flight constitute an N-body orbital mechanics problem. Figure 2.2 shows a typical formation system including N satellites orbiting a central body (Earth). The mass of satellite i is denoted by and that of the Earth is The satellites center of mass is located relative to the Earths center of mass by vector im eM thii R and relative to th j satellites center of mass by vector (rR ijr ijijR ). Individually, the satellites are modeled as spherically homogeneous bodies with the mass acting through their respective mass center (c.m.); however, the Earth is modeled as a distributed mass body. Thus, for each satellite i, the gravitational potential with respect to the Earth is given by Eq. (2.1) [51]. 221{1[()(cos)()(cossin)(cos)]}kkeekkkmkmkmkkmiiiRRUJPCmSmPRRR (2.1) Here i R is the magnitude of the satellite position vector i R is the equatorial radius of the Earth, and and eR are the latitude and longitude of the satellite. is the order Legendre polynomial function, and is the associated Legendre functions of degree kP thk kmP n

PAGE 33

15 with order The gravitational parameter m is 3.986, are the zonal harmonic coefficients, and and are the sectional harmonic coefficients for and the tesseral harmonic coefficients for 5310/kmsm 2kJ kmC kmS km k constant coefficients. Figure 2.2 N bodies system The first summation in this potential function is related to the latitude variation in the mass distribution of the Earth (see Figure 2.3a), and the second term is related to the longitude variation (see Figure 2.4b). The darker shaded regions represent areas of greater mass density than the lighter shaded regions. According to Prussing [52], due to the Earths spin, the longitudinal variation in the Earths mass distribution can be neglected for all satellites except geostationary satellites. Then the gravitational potential of any arbitrary satellite with respect to the Earth, is reduced to Eq. (2.2). Thus, for satellites in formation systems (which are non-geostationary), the tesseral and sectoral variations are neglected and the Earth is assumed axially symmetric; that is, the gravitational potential is i

PAGE 34

16 2[1()(cos)]kekkkiiRUJPRR (2.2) a Tesseral b Sectoral Figure 2.3 Zones of variations It is well known [52, 53] that the contributions of the zonal terms (J 2 J 3 ,) to the gravitational potential are on the order of 10 or smaller. Although they are small, they have significant influence on the precision of the orbit, which is a primary interest to formation flight. However, these contributions will be considered as perturbing acceleration and will be developed in Section 2.4. What follows is the classical Keperian orbit analysis. g3 ida 23JiiiiUuRR idRa (2.3) Neglecting the Earths zonal harmonic effects and considering inter-satellite attraction, the motion of the i satellite is governed by Eq. (2.4), where th 3jijijGmrr is the force exerted by th j satellite in the formation, and ida is the disturbing acceleration which consists of the zonal harmonic effects and other acceleration (e.g., from applied force).

PAGE 35

17 331331NjiiijidjiijjiNjiijidjiijjiGm R RrRrRraRr a (2.4) Typically, formations consist of micro-satellites whose masses m and relative distance are the on the order of 100 Kg and 10 Km. For a typical low Earth orbit scenario ( Km), j ijr 7000iR 2032226.672210/jjijijuGmkmsrr and 3220.0081/ikmsR ida For this typical LEO scenario, it is obvious that the contribution of the inter-satellite attractions are even smaller than those non-spherical mass distribution of the Earth, and so their contributions are also considered a disturbance and included in Thus, the N-body satellite formation problem can be simplified to an individual Earth-satellite 2-body relative motion problem as shown in Eq. (2.5). Here ida includes the normalized control efforts, disturbance and N-body attractions. 3iiidi R RaR (2.5) 2.3 General 2-Body Relative Motion Models In this section, the general 2-body dynamic models for the relative motion are derived for the formation flying system, which assumes a leader/follower or master/slave architecture. In the development that follows, each satellite is assumed to experience 2-body relative motion with respect to the Earth and there are no assumptions regarding to the eccentricity of the leaders orbit.

PAGE 36

18 2.3.1 Relative Motion Dynamic Model (Point Mass Case) Consider the 2-body relative motion model shown in Figure 2.4. L R F is the position vector from the Earths center of mass to the leaders center of mass, R is the position vector from the Earths center of mass to the followers center of mass, and (in Section 2.2, is used) is the relative position vector from the leader to the follower satellite. An Earth-Centered-Inertial (ECI) coordinate frame is attached to the Earth center as shown in The angles ijr , i and f are the right ascension, inclination, argument of periapsis, and true anomaly, respectively, of the satellite. Furthermore, it is customary to define as the longitude of periapsis in the singular cases, and u 0 f as the argument of latitude at epoch. Also, and h v are the specific angular momentum and satellites velocity, respectively. In Figure 2.6, local vertical local horizontal (LVLH) coordinate frame is attached to the leader as shown with x in the radial direction, is parallel to specific angular momentum (i.e., normal to the orbit plane), and is established by the right-hand rule. z y Figure 2. 5 Figure 2. 4 Relative dynamic model

PAGE 37

19 Figure 2. 5 Coordinate system of the relative position model. Figure 2. 6 Coordinate system of the leader satellites frame As defined in Figure 2.5, the position of the followers c.m. relative to the leaders c.m. is given by F L R R (2.6) As discussed above, each vehicles motion is defined by a 2-body model and is shown in Eqs. (2.7) and (2.8). Note that, here symbol a is replaced by u for the purpose

PAGE 38

20 of representing all kinds of external input accelerations, includes the disturbing acceleration due to un-modeled effects and the mass normalized control actions. 3FFFF R RuR (2.7) 3LLLL R RuR (2.8) Thus, the relation motion of the follower with respect to the leader is governed by 33[()](LFLLFFLLFR ) R RRRuRR u (2.9) Now, it is customary to coordinatize in the local frame of the leader. Thus, the second time derivative of with respect to the ECI frame expressed in the LVLH is 2(LLLLLLLLLLLL )L (2.10) where implies the derivative with respect to the ECI frame, and denotes the derivative with respect to a LVLH frame. The left superscript in the vector represents the frame in which the vector is coordinatized, and the right subscript denotes the quantity of interest (i.e., the leader or a follower). For example, and LL L L are the angular velocity and angular acceleration of the LVLH frame with respect to the ECI frame, expressed in the LVLH frame. Furthermore, L and L are the local relative velocity and acceleration expressed in the leaders frame. Substitution Eq.(2.10) into Eq. (2.9) yields the governing equation for the relative position of the follower with respect to the leader.

PAGE 39

21 332()[()](LLLLLLLLLLLLLLFLLLLFLFR )L R RuuRR (2.11) Equation (2.11) is the vector form of the generalized relative motion model for the point mass case. Now, to evaluate this expression, the quantities defined in the expression must be evaluated. For example, the angular accelerations of the leader and follower satellites, coordinatized in the LVLH frame, are 00TLLLf and where 00TLFFf L f and F f are the true anomaly accelerations of leader and follower, respectively. The position vector of the follower and leader satellites are defined to be 00 TFRR FF and 00TLLRR L To express F R in leaders frame, the transform matrix between the two frames must be determined. Using Figure 2.5, the position vector of the follower coordinatized in LVLH can be written as Eq. (2.12) and the rotation matrix from the followers coordinate system to the leaders coordinate system can be written as shown in Eq. (2.13), where i C represents the fundamental rotation matrix about the i axis as in Eqs. (2.14) and (2.15). th LFLF FF R CR (2.12) 00//313313()()()()()(TLLLFFFLFECILECIFCCCCuCiCCCiCu ) (2.13) 1100()0cossin0sincosCxxx x x (2.14)

PAGE 40

22 3cossin0()sincos000xxCxxx 1 (2.15) The remaining terms of Eq. (2.11) are TLxyz (2.16) 0yLLLxLL f f (2.17) 22()0xLLLLyLLL f f (2.18) 2220yLLLxLL f f (2.19) If the i element of the 1 th st column of L F C is denoted by b, then substitution of Eqs. (2.12) (2.19) into (2.11) yields the governing equation for relative motion. i 231323233332[()]2[()][()]LyLxLyLFxFxLxLLFLxLyLxFyFyLyLLFLFzFzLzLFR f ffRRbuRRR u f ffRbuRRRRbuuRR u (2.20)

PAGE 41

23 xFyFFzFuuuu (2.21) xLyLLzLuuuu (2.22) Eq. (2.20) is the generalized dynamics model for the relative motion in formation flight. It is well know that analytic time-dependant solutions of the classical 2-body problem cannot be obtained unless the orbit is circular. Thus, for the closed orbit case, the quantities L f and L f in Eq. (2.20) cannot be determined. It is also well known that the analytic solution resulting from the classical 2-body problem is dependant on an angular measurement, which is the true anomaly. Thus, in order to be able to evaluate L f and L f the governing equation expressed by Eq. (2.20) must be regularized. Eqs. (2.23) and (2.24) represents the leaders angular velocity and acceleration as functions of the leaders true anomaly L f semimajor and eccentricity La Le 222222(1)(1)(1cos)(1)LLLLLLLLLLaeaeeffrae 2 (2.23) 33232(1cos)sin(1)LLLLLLeefffae L (2.24) Using the chain rule to change the independent variable from time to the leaders true anomaly as the example listed in Eq. (2.25), the most generalized model for the relative motion in formation flight is achieved as given by Eq. (2.26) or in the vector form Eq. (2.27) [54, 55].

PAGE 42

24 '22''''2()LLLLLdfdgdggfdtdfdtdgdgfgfgfdtdt L (2.25) 2''''231322''''23232'''33321{2[()1{2[()1{[())]}LxxLyyLxLLFLLFLLyyLxxLyLFLLFLLzzLFLFLR ]}]} f fffRRbRRfRffffRbRRfRfRbRRf (2.26) ''''2331{2()[()LLLLLLLLLLLLLLLFLLLF ]}L f ffRRRRR (2.27) In Eq. (2.26), the terms [1,2,3]FiRbi can be calculated from the leaders position. This term is the coordinates of the followers position expressed in LVLH. Therefore, thix 1FLRbR 2Rb Fy and 3Rb Fz For the purpose of evaluating L f L f and the desired attitude as shown in Section 2.3.5, the followers six orbital elements: Tfiea i.e., semimajor, eccentricity, right ascension, inclination, argument of periapsis, and true anomaly (as shown in Figure 2.7) are need to be obtained. Therefore the relations between the position vector R absolute velocity vector and the six orbital elements are needed. v

PAGE 43

25 Figure 2. 7 Orbital elements In terms of the orbital elements, the specific angular momentum magnitude is )1(2eah the argument of latitude at epoch is uf 0 and the magnitude of the position vector is 2(1)1cosaeRef From these quantities, it is customary to write the position and absolute velocity as 00000coscossinsincossincoscossincossinsinECIuu i R Ruuui i (2.28) 00000[cos(sinsin)sin(coscos)cos][sin(sinsin)cos(coscos)cos](coscos)sinECIueueivueuehuei i (2.29) In order to evaluate Eqs. (2.28) and (2.29), knowledge of the orbital elements is required. They can be computed if the position and velocity vectors are known at some epoch condition. The six orbit elements can be found as described below. First, the semi-

PAGE 44

26 major axis can be determined from the position and velocity; the result is given in Eq. (2.30). Next, from the specific angular momentum, hR v the eccentricity vector and inclination are determined from Eq. (2.31) and Eq. (2.32), where is the coordinate of the angular momentum in the Earth frame. zh z semi-major axis 212(vaR ) (2.30) eccentricity 1[evh ]R R (2.31) inclination 1coszhih (2.32) The other three elements have singularity cases. Case 1: If the angles 0i , and f are not defined (see Figure 2.8). As such, only the true longitude Lf is calculated here. This case includes both and e. 0e 0 true longitude 11cos02cos0yyRxRRLRxRR (2.33) Case 2: If and the angles 0i 0e and f are not defined (see ). The nodal vector is calculated by n [0,Th 0,1]h n The right ascension and argument of latitude are obtained by Eqs. (2.34) and (2.35). Figure 2. 7

PAGE 45

27 right ascension 11cos02cos0xyxynnnnnn (2.34) true longitude at epoch 101cos02cosRnRvRnuRnRvRn 0 (2.35) Case 3 : when i and e 0 0 the right ascension argument of periapsis and true anomaly f are all defined and calculated from Eqs. (2.36) (2.39). nodal vector [0,0,1]Thnh (2.36) right ascension 11cos02cosxyxynnnnnn 0 (2.37) argument of periapsis 11cos02coszzneeneneene 0 (2.38) true anomaly 11cos02cosReRvRefReRvRe 0 (2.39) To here, the regularized dynamics model of the relative position for the formation system is derived. Considering the singularity problem, this model can be used under both circular and elliptic reference orbit cases. In the next Section, simulation results will show the correctness of the model.

PAGE 46

28 2.3.2 Simulation Results and Discussion (Point Mass Case) The following simulations are used to validate the generalized dynamic model developed in this dissertation. Responses from the regularized model are generated and compared with the Keplerian model solution as shown in Eq. (2.28). 2.3.2.1 Initial condition Time dependant model listed in Eq. (2.20): Initial relative position and relative velocity are calculated from Eqs. (2.40) and (2.41) below, where the subscript 0 represents the initial condition. Here 0FF R is the followers velocity with respect to followers frame expressed in followers frame, and LF denotes the follower frames angular velocity with respect to the leaders frame. 000000000//313313()()()()()(TLLLFFFLECILECIFCCCCuCiCCCiCu ) 000LFFFLCRR 0LL (2.40) 0000()()LFLFFLFFLFLCRCRR 00LFL (2.41) Regularized model listed in Eq. (2.27). Given a set of the initial condition for the six orbital elements Eqs.(2.28) and (2.29) are used to compute the initial absolute position and velocity 000000[,,,,,]Tifae 0 R and 0v Through the rotation matrices, the localized velocity and position with respect to time are calculated. The initial condition for relative velocity with respect to the leaders true anomaly is found from the chain rule in Eq. (2.25). 2.3.2.2 Results and discussion Due to the singularity cases discussed in the previous section, three scenarios are simulated. The first scenario involves a zero inclination and a non-zero eccentricity of the

PAGE 47

29 followers orbit, the second scenario involves a nonzero inclination and a zero eccentricity of the followers orbit, and the last scenario involves both nonzero eccentricity and inclination of the followers orbit without loss of generalization. In all the simulations, right ascension is selected to be zero, and the argument of periapsis is 90. o Table 2-1 shows all the cases that were tested. These cases include all combinations of the singularity scenarios. It should be noted that only the generalized relative position dynamics models are checked in these simulations. The dimension of the Earth has not been considered. Therefore, the satellite may go into the Earth, which is not practical. The derived generalized models are compared with the ideal Kepelerian model by the norm2 difference percentage || / K EPLERIANREGULARIZEDKEPLERIANRRR Table 21 Different cases been simulated. Case La (km) Fa (km) Li (deg) Fi (deg) L f (deg) F f (deg) Le Fe 1 6878 6878 0 5 270 265 0 0 2 6878 6878 0 5 270 265 0 0.5 3 6878 6878 0 5 270 265 0.5 0 4 6878 6878 0 5 270 265 0.5 0.75 5 6878 6878 0 5 265 270 0 0 6 6878 6878 0 5 265 270 0 0.5 7 6878 6878 0 5 265 270 0.5 0 8 6878 6878 0 5 265 270 0.5 0.75 9 6878 6878 5 0 270 265 0 0 10 6878 6878 5 0 270 265 0 0.5 11 6878 6878 5 0 270 265 0.5 0 12 6878 6878 5 0 270 265 0.5 0.75 13 6878 6878 5 0 265 270 0 0 14 6878 6878 5 0 265 270 0 0.5 15 6878 6878 5 0 265 270 0.5 0 16 6878 6878 5 0 265 270 0.5 0.75 17 6878 6898 0 5 270 265 0 0 18 6878 6898 0 5 270 265 0 0.5 19 6878 6898 0 5 270 265 0.5 0 20 6878 6898 0 5 270 265 0.5 0.75 21 6878 6898 0 5 265 270 0 0 22 6878 6898 0 5 265 270 0 0.5 23 6878 6898 0 5 265 270 0.5 0 24 6878 6898 0 5 265 270 0.5 0.75

PAGE 48

30 Table 2-1 Continued Case L a (km) Fa (km) L i (deg) Fi (deg) L f (deg) F f (deg) L e Fe 25 6878 6898 5 0 270 265 0 0 26 6878 6898 5 0 270 265 0 0.5 27 6878 6898 5 0 270 265 0.5 0 28 6878 6898 5 0 270 265 0.5 0.75 29 6878 6898 5 0 265 270 0 0 30 6878 6898 5 0 265 270 0 0.5 31 6878 6898 5 0 265 270 0.5 0 32 6878 6898 5 0 265 270 0.5 0.75 33 6898 6878 0 5 270 265 0 0 34 6898 6878 0 5 270 265 0 0.5 35 6898 6878 0 5 270 265 0.5 0 36 6898 6878 0 5 270 265 0.5 0.75 37 6898 6878 0 5 265 270 0 0 38 6898 6878 0 5 265 270 0 0.5 39 6898 6878 0 5 265 270 0.5 0 40 6898 6878 0 5 265 270 0.5 0.75 41 6898 6878 5 0 270 265 0 0 42 6898 6878 5 0 270 265 0 0.5 43 6898 6878 5 0 270 265 0.5 0 44 6898 6878 5 0 270 265 0.5 0.75 45 6898 6878 5 0 265 270 0 0 46 6898 6878 5 0 265 270 0 0.5 47 6898 6878 5 0 265 270 0.5 0 48 6898 6878 5 0 265 270 0.5 0.75 Simulation Result 1, (Case 1 in Tabel 2-1) : In this case, the eccentricities of the leader and follower orbit are zero, but the inclination of the follower is ( 5o 05LFff ), and all the other orbit information can be found in Table 2-1. Figure 2.8 depicits the orbits of the leader and follower satellites. In this figure, the follower satellite is in an inclined orbit and behind the leader. Figures 2.9 through Figure 2.11 show the relative position expressed in the LVLH frame. It is seen the direction period trend in Figure 2.8, where in half period, the follower is above the leader, and in the other half, below the leader. For the follower is behind the leader, in LVLH coordinates, 5o z 5o x and directions, the follower is always later and lower than the leader with periodic variations. Also the 3D relative position is shown in Figure 2.12. The 2-norm difference percentage shown in the results (Figure 2.13) is on the order y

PAGE 49

31 of that can be regarded as the MATLAB precision in numerical calculation. Therefore, for the circular orbit case, the regularized model gives results that are compatible with the Kepelerian model. 16(10)O 0 1x 104 1000 500 0 500 1000 y axis in E Figure 2. 8 Orbital information 0 2000 4000 6000 8000 10000 12000 -55 -50 -45 -40 -35 -30 -25 time second (2 period)x coordinates (km)Relative Position in x axis Figure 2. 9 Relative position in x direction 0 2000 4000 6000 8000 10000 12000 -610 -605 -600 -595 -590 -585 time second (2 period)y coor di na t es (k m ) Relative Position in y axis Figure 2. 10 Relative position in y direction 0 2000 4000 6000 8000 10000 12000 -600 -400 -200 0 200 400 600 time second (2 period) z c oor di n a t e s (k m ) Relative Position in z axis Figure 2. 11 Relative position in z direction -1 -0.5 0 0.5 1x 104 -1 x axis in ECI (km) Orbit Information CI (km) z axis in ECI (km)EARTH Follower Leader

PAGE 50

32 -60 -50 -40 -30 -20 -620 -600 -580 -1000 -500 0 500 1000 x axis in Local Frame (km) Relative position y axis in Local Frame (km) z axis in Local Frame (km) Figure 2. 1 Relative position in 3D 0 2000 4000 6000 8000 10000 0 1 2 3 4x 10-16 Time Seconds (sec)Norm 2 Difference PercentagesDifference: Kepelerian & Regularized Figure 2. 2 Model difference Simulation Result 2, (Case 9 in Table 2.1) : The eccentricity and inclination of the follower satellite are both zero. As shown in Figure 2.14, the regularized model has the difference on the order of The relative position results are not demonstrated here because they are similar to the previous case. 16(10)O 0 2000 4000 6000 8000 10000 0 1 2 3 4x 10-16 Time Seconds (sec)Norm 2 Difference PercentagesDifference: Kepelerian & Regularized Figure 2. 3 Case 9 model difference Simulation Result 3, (Case 20 in Table 2.1) :

PAGE 51

33 The eccentricity and inclination are all nonzero for the follower and leaders orbits, and the model difference is shown in Figure 2.15. The difference percentage is on the order of This relatively larger difference in comparison to the two previous cases is due to the additional computations required for th elliptic orbit of the follower. Unlike the circular case, the true anomaly of the follower in the elliptic orbit cannot be obtained from the leaders frame directly. In calculating the followers true anomaly for the Keplerian orbit, first the leaders true anomaly is determined from the simulation step. Next, the regularized model uses the leaders true anomaly to calculate the followers absolute position and velocity. Finally, the followers true anomaly is determined in Eqs. (2.42) and (2.43). Computations have been conducted, which shows that the difference is coming from those above numerical steps. Furthermore, the spikes in the plot are associated with the true anomalies of 0, 12(10)O 12(10)O and 2 etc. Therefore, the regularized dynamics model for the elliptic orbit case is also compatible with the Kepelerian model. FLRR (2.42) ()cos()ECIFxFFFF R fR (2.43) From the above cases, the regularized relative dynamics model is correct for the formation flying system. The 48 cases investigated in Table 2-1 indicate that the regularized model produce results are numerically exact to those of the time dependant model, and thus the model appears correct. In the next condition, the special cases of the relative dynamics model are investigated.

PAGE 52

34 0 2000 4000 6000 8000 10000 0 0.2 0.4 0.6 0.8 1 1.2 1.4x 10-12 Difference: Kepelerian & RegularizedNorm 2 Difference PercentagesTime Seconds (sec) Figure 2. 15 Case 20 model difference 2.3.3 Specialized Relative Motion Model (Point Mass Case) Generalized dynamic models of the relative position motion in the formation system have been developed in Section 2.3.1. In this Section, three specialized models are studied. The first two have been heavily used in applications and literatures which are the C-W equation and the in-plane nonlinear circular reference orbit equation. The third one is the in-plane nonlinear elliptic reference orbit equation. Here in-plane represents the case where the leader and follower satellites are in the same plane. 2.3.3.1 Linear model (Clohessy-Wiltshire equation) If the reference orbit is a circular orbit, the general 2-body relative motion dynamic model derived in Section 2.3.1 can be simplified to a linear system, which is the well-known Clohessy-Wiltshire (C-W) equation. According to Section 2.3.1, the second derivative of the relative position with respect to the ECI frame can be expressed in a LVLH frame, here, i.e., the leaders frame as Eq. (2.44). For the general 2-body motion 3FFFF R RuR let 3()gRRR FL then Eq. (2.45) can be achieved. Substitution of the relative position R R into Eq. (2.45) and using a Taylor series expansion about

PAGE 53

35 reference position L R ) Eq. (2.46) can be derived. Neglecting the higher order terms and utilizing (LLL R gRu the linear model is derived as Eq. (2.47). Now the equation is written in the reference frame, i.e., leaders frame, using Eq. (2.48) and FLFLuuu)L LL FLu 000TLL 00LLf T L 001 2(LLLLLLLLLLLL (2.44) ()FFF R gRu (2.45) 2()()()LLLFLgR R gRuR (2.46) ()LFLgRuuR L (2.47) ()2()LLLLLLLLLLLLgRR (2.48) In the circular reference frame case, , 00LTf L 00TLRR and 3200100LRR ()LLgR Finally, the C-W equation is derived as Eq. (2.49), where matrices and A B are shown in Eqs. (2.50) and (2.51).

PAGE 54

36 (2.49) xxxyyyzzzuABuu FLxFLyFLz 2230000000LLfA f (2.50) 0202000LLfBf 00 (2.51) The C-W equation was derived in the nineteenth century [52], which originally described the motion of the Moon relative to the Earth. Although it is different from the formation flying application, they both consider small displacements relative to a known reference motion. This linear model has been used extensively in formation flight models [19, 20, 21, 22, 23, 24, 26, 27, 28, 29, 30]. 2.3.3.2 Circular reference orbit nonlinear model Due to the uncoupling between the in-plane and the out-of-plane motion, nonlinear co-planar relative position equations with respect to a circular reference orbit will be shown in this part. In this case, the leader satellite is limited in the circular orbit. Therefore, proper coordinates are set to let the inclinations for the leader and the follower satellites to be the same, i.e., ii 0LFi Here the right ascension is assumed to be zero. Based on these conditions, from Eq. (2.20), the circular reference orbit in-plane model can be simplified to Eq. (2.52), which has the same results as derived from the Keplerian model.

PAGE 55

37 230000323000032[()(coscossinsin)]2[()(sincoscossin)]LxLyLFLFLFxFxLxLLFLyLxFLFLFyFyLyLLFR f fRRuuuuuRRR u f fRuuuuuRR u (2.52) 2.3.3.3 Elliptic reference orbit nonlinear model Just as shown in section 2.3.1, in elliptic reference orbit case, it is not convenient to represent the model as a function of time, because the angular velocity of the leader satellite is not a constant. Therefore, the true anomaly is used as the independent variable. Based on the dynamic model derived in the previous section and using the chain rule dtdfdfdgdtdg the governing equations for the in-plane elliptic reference orbit relative motion can be re-written as Eq. (2.53), where the notations are listed in Eqs. (2.54) (2.57). As Section 2.3.3.2, this model is also an in-plane special case derived from the general relative dynamics model. ,,,2,213/2212,,,2,223/22121{2()[()]()()1{2()[()]()()xLxLyLxLxLxxFLyLyLxLyLyLyyFL }}xLyL f fffcRuddf u f fffcRuddf u (2.53) 222212[(1)](1cos)LLLxLyLLaedRRef (2.54) 2222(1)2(1cos)LLxyLLaedef x (2.55)

PAGE 56

38 23(1)LLdhae L (2.56) 3/23/2121213/23/2112[()]()ccdddddd (2.57) Up to now, three special cases of the relative motion models have been introduced. In the next section, how precisely the C-W equation is are shown by simulation. 2.3.4 Simulation Results and Discussion (Point Mass Case) Simulation results show the difference between C-W equations and the generalized relative position model. Results demonstrate the necessary of using nonlinear model instead of the C-W linear model in the formation flight system. The in-plane circular and elliptic reference orbit model are not compared, since it is identical to the generalized relative position model under the in-plane assumption. The leader is in an equatorial orbit with 6878 km radius and 0 true longitude. The follower is also in the equatorial orbit, with 0 true longitude. But the followers radius is 6888 km. Initial relative velocity with respect to the leaders frame is calculated by o o 00flFLvv R R Figure 2.16 and Figure 2.17 show the relative position results from the C-W equation and the generalized relative dynamics model. In Figure 2.18, the difference between two models is shown. The difference is calculated by the 2-norm between these two models at each time second. C-W equation has analytic solution as shown in Eq. (2.58), where Lnf 0i and 0i are the initial relative position and velocity for Although C-W equation is easier to calculate, it can only be used in the [,,]ixyz

PAGE 57

39 circular reference orbit. Furthermore, due to the truncation of the high order terms, the solution is not periodic, which can be seen from Eq. (2.58). These are the reasons why the nonlinear model is adopted in this dissertation. 000000cossin2(43cos)sin(1cos)24sin6(sin)(1cos)xxxxyxyzxyxntntnntntntnnntntntntntnn 003y (2.58) 0 0.5 1 1.5 2x 104 -60 -40 -20 0 20 40 60 Time Seconds ( 3 Periods ) Relative Position in xComparison of CW and Generalized Model CW Generalized Model Figure 2. 16 Relative position in x direction 0 0.5 1 1.5 2x 104 -900 -800 -700 -600 -500 -400 -300 -200 -100 0 100 Time Seconds ( 3 Periods ) Relative Position in yComparison of CW and Generalized Model CW Generalized Model Figure 2. 17 Relative position in y direction 0 0.5 1 1.5 2x 104 0 20 40 60 80 100 120 Time Seconds ( 3 Periods ) Difference (km)Difference Between CW and Generalized Model Figure 2. 18 Difference between CW and generalized model

PAGE 58

40 2.3.5 Relative Motion Dynamic Model (Rigid Body Case) In section 2.3.1, the satellites are regarded as point mass vehicles and the relative position model is derived. In this section, the rigid body case is studied, which includes both position (orbit) and orientation (attitude) motions. To date, most literatures [19, 20, 21, 22, 47, 48, 49] regard the orbit and orientation of the formation flight as an uncoupled system. The relative position and attitude motions are studied independently. However, this is not necessary true. For example, missions as EO-1 and LandSat7, which were required to collectively point to a single spot on the surface of the Earth, once or twice a day. Furthermore, one of the final objectives of the formation flight system is pointing and positioning managed collectively. Therefore, the coupled position and attitude models must be studied. In this dissertation, the objective of pointing managed collectively is simplified to be the following. All the satellites in the formation flight system are required to point to the nadir of the leader satellite. Also, all the satellites are assumed to be gravity gradient stabilized and thus are pointing to their own nadir points under the passively stabilized condition. Based on the above simplification and assumption, only the desired orientations of the follower satellites are coupled with the real time relative position. In this section, first, the desired attitude model is derived. Then, the attitude governing equation is studied. For the convenience of using sliding mode controllers, which will be discussed in Chapter 3, the Euler angle representation for attitude are adopted. 2.3.5.1 Desired attitude For a gravity gradient stabilized satellite, the body frame is chosen as shown in Figure 2.19. For instance, shown in the followers body frame, axis points to the z

PAGE 59

41 followers nadir point FN x is perpendicular to axis z and in the followers orbit plane, and is determined by the right hand rule. The desired attitude of this follower satellite is that the direction points to the nadir of the leader satellite y z L N Define two error angles, an elevation error EL along latitude and an azimuth error AZ along longitude. In order to calculate these error angles from the real-time relative position, three theories coming from the spherical geometry are used. B C A B RR Figure 2. 19 Satellite attitude body frame Theorem 1 [Spherical Pythagorean Theorem Ref. 56] In a right triangle on a sphere of radius R with a right angle at vertex C as shown in Figure 2.20, angles ABC and are facing sides AC and AB separately. The lengths of sides B C AC and are a, b and c, then cos()cos()cos()ca b R (2.59)

PAGE 60

42 Theorem 2 [Spherical Sine Theorem, Ref. 56] Let ABC be a spherical triangle on a sphere of radius R. Let a, b, and c denote the lengths of the sides, and let, A B and denote the interior angles at each vertex. Then C sin()sin()sin()sin()sin()sin()ab c R RAB RC (2.60) Figure 2. 20 Spherical right triangle Theorem 3 [Ref. 56] Suppose that A BC is a right triangle on the sphere of radius R with a right angle at C. Then coscos(/)a sinAB R Figures 2.21 and 2.22 show the definitions of the error angles. Let be the inclination difference between the follower and leader orbits. The curve line Fiii L FN A is on the longitude line of the follower satellite and elevation error angle F L is facing this arc line. Point A is the intersection of FNA and the leaders orbit plane. Along the latitude line of the leaders orbit, is the nadir point of the leader and azimuth error angle faces the arc line In Figure 2.21, the angle between ON LN LAN F and is OAF L In Figure 2.22, the angles between OB and OA OB and ON L are defined to be F and L separately. In order to obtain the desired attitude, the follower satellite first rotates

PAGE 61

43 along its longitude from the point to point FN A about its body axis x through Then, points to through angle FL z LN L F after rotating about its body axis In the the angles definitions are summarized below. These angles are to be used in calculating the error angles. y Figure 2. 23 Figure 2. 21 Desired attitude I

PAGE 62

44 Figure 2. 22 Desired attitude II Angles Definition Restriction 1F L Angle between lines and ON OD F 1FFiLi F 2FL Angle between lines OD and OA 2LFiLi L FL Angle between lines and FON OA 12FFFLLL u 0F Angle between lines OB and ON F 002Fu F Angle between lines and OB OA 02F 1F Angle between lines and OB OD 02 1F Figure 2. 23 Desired attitude III Now, the error angles are to be calculated from the real time relative position information. When the relative position of a follower satellite (i.e., x y and z ) is found, the orbit elements of this follower, such as and 0Fu Fi F can be determined as discussed in Section 2.3.1. In Figure 2.20, using the spherical sine theorem, in triangle

PAGE 63

45 F B DN can be found through Eq. (2.61). Also, in this triangle, 1FL 1F can be calculated by spherical Pythagorean theorem through Eq. (2.62). Note that, the quadrant of the cosine function is determined by the true anomaly u. 0F)Fi1F F F 2 2F F L 2)F (2.61) 110sin(sinsinFLu F0F 11010110cos(cos/cos)02cos(cos/cos)2oFFFFFifuLuifuLu (2.62) Here, needs to be added into F 1F because the follower orbits right ascension of the ascending node may change due to the control effort and the orbits of the leader and follower may not cross at the same point. Therefore, let 21FF Based on the theorem 3 and spherical sine theorem, in spherical triangle D BA and FL F is given by Eqs (2.63) and (2.64). Therefore, is obtained by F L 1FFLLL From these two angles and F the error angles are to be derived. (2.63) 12tan(sintanFL Li 102210220cos(coscos)22cos(coscos)oFFFFFFFFifuLifuL (2.64) Two rotations should be taken to achieve the desired attitude. According to the Figure 2.24, first, follower satellite rotates EL about the body axis x and second it rotates AZ about the body axis The following steps are used to generate the real time EL and AZ. All vectors are expressed in ECI frame for the convenience of calculation. y

PAGE 64

46 Figure 2. 24 Error angles in the coordinate system. EL is the elevator error angle, whereas, AZ is the azimuth error angle. ()313()()()00eECILLF R OARRiR (2.65) ()313()()() L xECILLLyzRFORRiR (2.66) ()()()(ECIECIECIFAFOOA ) (2.67) ()()1()()cos[]|||ECIECIECIECIFAFOELFAFO | (2.68) Considering about the rotation direction when the follower satellite is in different quadrant, the satellite rotates about x axis in negative direction when the true longitude of the follower is in [0 ,) and positive when [ ,2) Now, the error angles between actual and desired attitude have been derived as shown in Eqs. (2.69) and (2.72).

PAGE 65

47 ()()1()()()()1()()cos[]2||||0cos[]||||ECIECIECIECIFECIECIFECIECIFAFOELifFAFOFAFOifELFAFO (2.69) ()313()()()00eECILLLL R ONRRiR (2.70) ()()()ECIECIECILFNFOON L (2.71) ()()1()()cos[]|||ECIECILECIECILFAFNAZFAFN | (2.72) After the error angles between the desired attitude and actual open loop attitude are calculated, the desired attitude can be obtained through the direction cosine matrix form. For the gravitational gradient stabilized satellite, the unforced attitude only has one principal axis rotation, which is the rotation against the body axis with angular velocity, which is related to the orbit period. Lets denote this open loop Euler angle to be y op Then the desired attitude direction cosine matrix is calculated in Eq. (2.73), where subscript i means the rotation with respect to body axis. [1,2,3] thi 212()()(opDesCCAZCELC ) (2.73) 2.3.5.2 Attitude model The attitude motion of a rigid body is governed by BODYBODYBODYBODYJJ T (2.74) where J is the centroidal inertial dyadic of the body, BODY is the angular velocity of the body coordinatized in the body and T BODY is the extErnally applied torque

PAGE 66

48 (coordinatized in the body). In terms of Euler angles, the angular velocity can be written as Eq. (2.75), where the rotation squence is chosen to be 1-2-3 through 312[,,][,,]TTELAZ 3 Rotation 1 is for the elevation error angle and rotation 2 is for the azimuth error angle. Rotation 3 is used here only a supplement in order to get three differential equations and 3330 23cossinsinco000 The angular acceleration is derived as shown in Eq. (2.76). 11223333333323233cossin0sincos000sincosBODYxyzz 1212232320cos001sincoscosscossincsin 33332312323123312cossinsincos00sincoscoscoscossinsin ()S 2323coscossincossincossin0(BODYS 3123223001 1)(,)S S 332232200cos0sins0010011sin0cos 31323in0os001 123()S (2.75) 10 3312231323231312231323122cossincoscossinsinsinsincoscoscos (2.76) where is a trigonometric function of the last two Euler relations and is a column matrix of the Euler rates. Substituting Eq. (2.75) and (2.76) into Eq. (2.74) results in

PAGE 67

49 11()[BODYJSJSSJST ] (2.77) The attitude model in Eq. (2.77) is used for the large angle case. In the small angle perturbation, the infinitesimal angles can be used to linearize the attitude governing equations. For any small disturbance 123T the actual angular velocity and acceleration of the follower body can be written as 0BODYBODY and 0BODYBODY where subscript 0 represents the stabilized point. Neglecting the high order term, the small angle perturbation equations can be written as the following linear system CD T Here, C and D are constant matrices determined by the stabilization point and T is the control torque. Furthermore, because the relative position model is in terms of true anomaly instead of time, the attitude governing equation also needs to use the chain rule and be changes from the time dependant version to the true anomaly dependant version. 2.3.6 Simulation Results and Discussion (Rigid Body Case) Simulation Result 1: The leader satellite is in an equatorial circular orbit, with a 500 km altitude and the true longitude to be The follower satellite is also in a 500 km circular orbit with a inclination and a true longitude initially. Assume the satellites are gravity gradient stabilized and thus naturally points to their respective nadir. The open loop simulation results are show in Figure 2.25. Figure 2.25 shows that the elevator error angle is periodic. This is expected as shown by the angle definition in Figure 2.24 and the chosen orbits. The follower satellite rotates about 0o5o 5o x in the negative direction in the first half of the period, and then positive in the second half of the period. The azimuth error

PAGE 68

50 angle is also periodic. The maximum value in the elevator error angle can be shown to be on the order of 0.87 radians and this is confirmed in Figure 2.25. Simulation Result 2: In this scenario, the case when the inclination of the leader is non-zero is evaluated. The only differences in initial condition from the simulation 1 are the inclination of the leader is 5 and that of the follower is 10. The simulation results are shown in Figure 2.26 and are identical to those obtained in simulation case 1. Because the relative inclination between the leader and the follower has not been changed, thus the results are not changed. o o 0 2000 4000 6000 8000 10000 12000 -1 -0.5 0 0.5 1 Simulation Results 1Theta1 and Theta2Time Seconds ( 2 Periods ) Theta1Theta2 Figure 2. 25 Theta 1 (EL) and Theta2 (AZ) in radians 0 2000 4000 6000 8000 10000 12000 -1 -0.5 0 0.5 1 S i mu l at i on Resu l ts 2Theta1 and Theta2Time Seconds ( 2 Periods ) Theta1Theta2 Figure 2. 26 Theta 1 (EL) and Theta2 (AZ) in radians Simulation Result 3: The case when the follower is in an elliptic orbit is shown in Figure 2.28. In the point mass case, the Earth dimension ha not been considered. Here in the rigid body case, the Earth dimension is included. So, the elliptic orbit needs to be chosen carefully to ensure that the satellite will not fly into the Earth. The semi-major axis of the leader and the follower is chosen to be 13868 Km. Inclination of the leader is 0, and that of the

PAGE 69

51 follower is 5. Also the eccentricity of the follower satellite is 0.5. True longitudes of the leader and follower are and 5 respectively. As shown in Figure 2.28, o 0o o 1 goes to a large angle when the leaders true longitude arrives k where In the elliptic case, the error angles are also periodic. 0,k 1,2... x 1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 Ti ) a1a2 kg 0 1 2 3 404 Simulation Results 3me Seconds ( 2 Periods Theta1 and Theta2 Figure 2. 27 Theta1 (EL) and Theta2 (AZ) ThetThet 2.4 Perturbation Until now, all the equations and simulations have not included disturbance effects. Actually, a satellite in an Earth orbit experiences perturbations, although the absolute value is small. As stated in Prussing [52], The perturbations come from the lack of spherical symmetry of the Earth, the attraction of the Moon, and Sun, and atmosphere drag if the satellite is in the low orbit. The perturbation forces are normalized and expressed as accelerations (force per unit mass) and coordinatized in LVLH frame. In this dissertation, only perturbations associated with the atmospheric drag and the effects of the Earth oblateness are considered. KN/

PAGE 70

52 2.4.1 Atmospheric Drag The drag force acts in a direction opposite to the velocity of the satellite. For a satellite is in a low Earth orbit (LEO), it experiences a drag due the Earth atmosphere which can be represented as Eq. (2.78) [52]. 212DadCAam vv (2.78) Here, is the atmosphere density, C is the drag coefficient, D A is the cross-section area of the satellite, v and v are the velocity magnitude and unit vector of the satellite, and is the mass of the satellite. The quantity C is called the ballistic coefficient of the satellite. Subscript ad represents atmospheric drag. Therefore, the atmospheric drags for the leader and the follower are m mAD/ 212DLLLLCAav adL v m and 212DFFadFFa FCAm v v respectively. The relative atmospheric drag for the generalized motion in a formation system is calculated by Eq. (2.79), where ti 1,2, 3 i is determined by LFCadFaadL , and a. Here, are the leaders LVLH unit vectors. The table below listed the coefficients for the atmospheric drag model. The altitude of the satellite is on the order of 500 Km in the simulation. The dimension and weight of the satellite are chosen from ST-5 [51]. Lree Le ]TzLe [,,L Atmospheric density Drag Coefficient D C Area A Mass m 1334.8910/kgm (500km altitude) 2 (typical) 20.2860.542m 19.5kg 123rLLzLadFadLaatetete (2.79)

PAGE 71

53 Figure 2. 28 Atmosphere drag 2.4.2 Earth Oblateness As was discussed in Section 2.2, the averaging effects of the Earths spin results in a gravitational potential of the satellite with respect to the Earth as expressed in Eq. (2.80) [52]. R is the distance from the Earth center to the satellite, e R is the Earth equatorial radius, is the Legendre polynomial function of the kth order and is the colatitude. The first three coefficients are kP 00108263.02 J , and. 00254.03J 000 00000161.04J 2[1()(cos)]kekkkRUJPRR (2.80) Since is significantly bigger than the others, the satellites potential with respect to the Earth is approximately as 2J

PAGE 72

54 222[1()(cos)]eRUJPRR (2.81) where the Legendra polynomial is 221(cos)(3cos1)2P Can be determined as Figure 2.29, the quantity cos is simply Z R where Z and R are Z component and the magnitude of the satellites position respectively. Thus, the potential including oblateness effects is 22221{1()[31)]}2eRZUJRRR (2.82) In the expression for the potential, the first term is related to the central force, which has been previously considered. [52] Thus, second term is the disturbance related to the Earths oblateness 2222321[(31)]2eJRZUJRR (2.83) Figure 2. 29 The perturbation stems from the Earth oblateness

PAGE 73

55 The perturbating acceleration can be derived as Eq. (2.84) [52], where the subscript represents the Earths oblateness. In order to accommodate the perturbations into the general dynamics model for the relative position in the formation system, the difference between the follower and the leaders oblateness pertubation need to be derived as shown in Eqs. (2.85) and (2.86). Therefore, the term perturbation for the relative position model is written as Eq. (2.87) where obl 2J1, 2,3ipi is determined by Eqs. (2.85) and (2.86). 22222224653153()22JJroblerUUaeZRZ e Z ZJReJRZRRR (2.84) 2224651533()22FerFeoblFFFFZaJReJRRRR 22FZZ (2.85) 2224651533()22LerLeoblLLLLZaJReJRRRR 22LZZ (2.86) 123rLLzLoblFoblLaapepepe (2.87) Include the and the atmospheric drag terms derived in Section (2.4.1), the generalized relative position model derived in Section 2.3 becomes Eq. (2.88). 2J

PAGE 74

56 22''''31113222''''322232'''3333321{2[()]1{2[()]1{[()]}LxxyyxLFLLLLLFLLyyxxyFLLLLLFLLzzFLLFLR }} f fffRRbtRRfR p f fffRbtpRRfRfRbtpRRf (2.88) 2.4.3 Perturbation Effects on the Generalized Dynamics Model Simulation 1: In this section, the effects of the perturbations are shown. The satellites orbit is a circular orbit with 0-degree inclination and 500 km altitude. As shown in Figure 2.31, the altitude of the satellite is decreasing due to the atmospheric drag. Seen in Eq. (2.78), the atmospheric drag always acts on the velocity direction and decrease the velocity due to the negative sign. From the inverse law model (or Keplerian model), less velocity causes lower orbit and degrades the altitude of the orbit. However, the J 2 effect causes periodical variations on the satellite altitude. J 2 comes from the Earth oblateness, which is caused by the non-uniform distribution of the Earth density. But the total mass does not change. Therefore, satellites altitude changes periodically, but the maximum altitude maintained (as shown in Figure 2.32). For the simulation condition, where the J 0i 2 acceleration obl a can be simplified as 22432exJReR which is in the negative radial direction. This acceleration causes the satellite to be driven inside. Furthermore, J 2 changes right ascension and argument of perigee as shown in Eqs. (2.89) and (2.90). Therefore, the right ascension node and argument of perigee change due to the J 2 effects. In order to easily show the effects using simulation, J 2 coefficient is

PAGE 75

57 made to be 100 times larger than the correct number. Figure 2.33 shows the combined effects from the above two reasons. 2147/2222.0647410(cos)(1)Jaie (2.89) 2147/22221.0323710(45sin)(1)Jai e (2.90) Figure 2. 30 Atmospheric drag effect Figure 2. 31 J2 perturbation effect Figure 2. 32 Combined perturbation effect

PAGE 76

58 Simulation 2: In this part the effects of disturbance on the relative dynamics (Eq. (2.26)) model are shown in Figures 2.33 and 2.34. The leader satellite is in a 6878 km equatorial orbit. The follower satellite is the same radius circular orbit as the leader, except with 5 inclination. The true longitude of the leader and follower are and 5, respectively. The 2-norm difference between the case without perturbations and the case with perturbations are shown. The absolute errors are very small which are on the order of in the relative position model and in the relative attitude model. Since the errors are increasing with the time, these perturbations must be considered. o 0o o 4(10) 6(10) 0 2000 4000 6000 8000 10000 12000 0 1 2 3 4 5x 10-4 Perturbation's Effects on the Relative PositionTime secondError (km) Figure 2. 33 Error in position 0 2000 4000 6000 8000 10000 12000 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8x 10-6 Perturbation's Effects on the Desired QuaternionTime secondError Figure 2. 34 Error in attitude 2.4.4 Perturbation Effects on Keplerian Elements In the previous section, simulation results have shown the perturbation effects on the absolute and relative dynamic models. In order to access the perturbations effects on the six Keplerian orbital elements, the Gaussian Form of Lagranges Planetary Equations

PAGE 77

59 is derived according to reference [57]. Here, the brief derivation and explanation are illustrated. The orbital elements vector is listed in Eq. (2.91) and the detailed information is given in Table 2-2. [TiaeM ] (2.91) Table 22 Classical elements for an Earth satellite Quantity Specified by Symbol Definition Special Cases Reference Time Epoch 0T Time for which elements are specified Size Semimajor axis a Half the long axis of the ellipse Shape Eccentricity e Distance for center of ellipse to focus divided by semi-major axis; dimensionless ratio. 0: circular orbit Inclination i Angle between orbit plane and earths equatorial plane 0o :equatorial orbit; 90o : polar orbit Orbit plane Right ascension of the ascending node Angle measured at the center of the earth in the equatorial plane from the vernal equinox eastward to the ascending node. 0o : the ascending node is at the vernal equinox. Orientation of orbit in the plane Argument of perigree Angle measured at the center of the earth in the orbit plane form the ascending node to perigee. Location of the satellite in its orbit Mean anomaly at epoach M 360/otPt where is the period and P is the time difference between the epoch and the last perigee passage before epoch 0, when the epoch time is perigee. In order to apply the Lagranges Planetary Equations (LPE), the mean anomaly M must be expressed as M n where 3na is the mean motion, and is the time of pericenter passage. The rotation matrix from the Earth frame to the local (i.e., the satellite) frame is given by

PAGE 78

60 123123123coscossinsincoscossinsincoscossinsinsincoscossincossincoscoscoscoscossinsinsincossincoslllCmmmnnniiiiii iii (2.92) Then the satellites absolute position and velocity vectors can be written as Eq. (2.93) and Eq. (2.94). cossin0aEeRCbE (2.93) sin/(1cos)cos(1cos)0afEeEvCbfEeE (2.94) where E is calculated by 1ntan21Eee ta. 2f Assume the disturbance scale function is W. The system equations of the 2-body dynamic motions under perturbation can be written as Eqs. (2.95) and (2.96). Let the homogenous solution for the system equations to be Eqs. (2.97) and (2.98). Then the basic method of variation of parameters is used here, the solution of R and v are assumed to have the forms shown in Eqs. (2.99) and (2.100).

PAGE 79

61 dRvdt (2.95) 3[]TdvWRdtRR (2.96) (,) R Rt (2.97) (,)vvt (2.98) (,()) R Rtt (2.99) (,())vvtt (2.100) Taking the derivative of the vector R and v and substituting these derivative into Eqs. (2.95) and (2.96), the partial differential equations are derived as shown in Eqs. (2.103) and (2.104). dRRRddttdt (2.101) dvvvddttdt (2.102) 0Rddt (2.103) []TvdWdtR (2.104) Using []Tv left multiplied with (2.103) subtracting [ ]T R left multiplied with (2.104), the LPE are derived as Eq. (2.105)

PAGE 80

62 []TdWLdt (2.105) where [][]TWvvL TR In the mean time, reference [57] gives the derivative of the position and velocity vectors with respect to the six orbital elements as Eq. (2.106) and (2.107). 13133212121121212(1)(1)sin03(1)(1)2/(1)mlRRaelaeminllRRbntaememmaqnnllRRambemenn 222ln (2.106) 21211231312312222/(1)/(1)03cos/(1)2(1)2(1)(1)mlvvfbelfbemnllvvft 222lbf f bemmiaennlvfavfmeben maen 1121(1)lamen (2.107) Combined the LPE with above partial derivative equations, the governing equations are derived as Eq. (2.108). Notice that the matrix L is nonsingular if e or and i. Also, according to reference [57], the disturbance scalar functions partial derivations with respect to the six orbital elements are listed in Eq. (2.109). 0 1e 0

PAGE 81

63 32234411sinsinsincos2sin2dWdiWidtidtfabifabifabidiWbWdaWdtiedtnafabifaedebWbWdWbWdtdtaefaefaefafae cosW (2.108) 223sin203cos(1cos)2cossin0cos00(1)ssin000TTddTTddTTddRaftefabWWaftuRiuefabRiafWWaRuuaiebraeWWuru inf 2sin(1cos)0fbaefb (2.109) Finally, Gauss form of LPE is derived to be Eq. (2.110), where )1(2eap 1cospRef and ph From these system equations, it is easily to find that the perturbation force Tddzuu drd uu affects the six orbital elements. Therefore, the orbital elements are not constant under the perturbations.

PAGE 82

64 22sin()00sincos()00cos()sincossin()sin2sin20sin()cos0(cos2)()sin0drdRfhiRfhipfpRfRifudhehehiuadtaefaphhRepfpRfreMhhpfRebbpRfaheahe00000dzu f (2.110) 2.5 Conclusion The generalized dynamics model has been developed with the inclusion of perturbation effects. The developed dynamic model has been shown to reproduce the special case models currently in the literature. The derived generalized dynamics model can be used in any closed reference orbit (i.e., circular or elliptic reference cases). It is shown that for the scenarios considered in this dissertation, the desired attitude is coupled with the real time relative position. Perturbation effects from the atmospheric and Earths oblateness are studied and their effects on the Keplerian model, relative position, and orbital elements are evaluated. Although those effects are small in the magnitude sense, their cumulative influence is significant and needs to be considered.

PAGE 83

CHAPTER 3 CONTROL STRATEGIES: ORBIT 3.1 Introduction In Chapter 2, the generalized dynamic models, including both the point mass case and the rigid body case, have been intensively investigated. In Chapters 3, 4, and 5, the controllers are designed in order to maintain and/or re-configure the formation. The controllers for the formation system can be treated as a two-level architecture, in which the higher level dealing with the systematic control and the lower level related to the individual leader-follower pair. The higher level of the controllers will be studied in Chapter 5 where the information about how to manage the whole formation system is discussed. In Chapters 3 and 4, the lower level of the controllers is discussed for the leader and follower satellites pair to maintain and track the desired relative position and attitude. First, in this Chapter, controllers are designed for orbit maintenance; then in Chapter 4, attitude controllers and coupled attitude and orbit controllers are designed. In the literature, numerous controllers have been proposed for the relative position and attitude for the leader-follower formation flying system [21, 22, 26, 31, 40, 41, 47, 49, 58, 64, 65, 66]. Roughly, these controllers can be classified into either exact-model-based or uncertain-model-based controllers. In the exact-model-based controllers, Linear Quadratic (LQ) and Proportional-Integration-Differential (PID) are often used. It appears that the LQ controller is the first controller to have been applied to the formation flying system [21, 22, 26, 40, 58]. Based on the Clohessy-Wiltshire (C-W) model, Vassar and Adams [21, 22] proposed the 65

PAGE 84

66 continuous LQ controller in order to keep the relative position within the formation system. In order to design a controller for this relative position model that uses a pulse based thruster, Kapila [26] made some variations to the C-W model. One was discretizing the continuous model, and the other is using v instead of normalized forces. Through these changes, the C-W equation was developed as 1ii i x x v where i and 1i represent discrete time steps. and are constant matrices. Based on this discrete model, the pulse based controller was designed. Other exact-model-based controllers found in the literature are the PID class of controllers. For example, based on the C-W model, Robertson [40, 58] proposed a PD controller augmented with a thruster mapper. The controller has a format as shown in Eq. (0.1), where is the thruster mapper and is defined as shown in Eq. (0.2). and are the proportional and differential gain and ()gy pK dK x is the state error. is a performance criterion value and is a constant thruster force level. Additionally, Robertson [58] proposed a bang-off-bang controller, which is also based on the C-W model and has a similar idea to the one presented in reference [40]. y thrusterf ()()pdugyKxKx (3.1) (3.2) 21.510.75()10.7520.150thrusterthrusterthrusterthrusteryfyfgyyfyfotherwise All of the above exact-model-based controllers are based on the linear C-W model, which has some limitations in representing the true motion. They all assumed that the system model is well defined and all the states are known. Since, this assumption is

PAGE 85

67 not always realizable, and the controllers based on this assumption are sensitive to the spacecrafts parameter uncertainties, unexpected disturbance, and the initial attitude/position rates, the resulting performance will be less then desirable. Therefore, uncertain-model-based controllers [31, 41, 47, 49, 64] have appeared in the literature recently. For example, adaptive, Lyapunov-type and sliding mode controllers are now widely used for the system with model uncertainty. Bernstein [47] proposed an adaptive asymptotic tracking strategy for the spacecraft attitude of the formation system in the case that the spacecrafts inertia matrix is not known. By re-writing the attitude governing equation, the large-angle-tracking problem becomes an attitude-keeping problem near the desired point. Using this proposed adaptive law, the tracking error can be controlled to zero. Whats new in Bernsteins contribution is that, at the same time, the spacecrafts inertia matrix can be predicted. However, in Bernsteins work the attitude and position are assumed to be uncoupled, whereas in this Chapter, the controller will be designed for the coupled system. A Lyapunov-type controller is another typical method used for the uncertainties model, which is mainly applied to the relative position control problem. For instance, based on the circular reference orbit, Kapila [31] modeled the uncertainty as where and represent the mass and the normalized acceleration of a satellite separately. [,,/,,,]TFFeFLdxdydzmmMGmmaaa m a e M and are the Earth mass and the gravitational constant. Subscripts and are the leader and follower, whereas, G L F x , and represent the coordinates and disturbance. Then the model error is y zd where is the predicted model parameters. Accommodating this uncertainty, the

PAGE 86

68 Lyapunov-type controller is designed. This controller guarantees global asymptotic convergence of the position tracking error in th e presence of unknown, slow varying satellite masses, disturbance forces, and gravit y forces. Nevertheless, this method need to be extended for the nonlinear coupled system. Kim and Choi et al [49, 64] used the s liding model control method in their work on spacecraft attitude and the smart structure co ntrol. Using this sliding mode controller, the uncertainties related to the natural fre quency and external disturbances of the materials can be rejected. The sliding model w ill be demonstrated in this Chapter to have the capability of controlling the formation system well. All of the above uncertainty-model-based controller examples demonstrated the capability of system control under parametric uncertainties, such as mass, inertia, and gravitational constant. However, there are some other uncertainties and limitations that have not been considered; for example, un-modeled dynamics, actuator saturation and actuator resolution. Furthermor e, the full nonlinear model should be used instead of the linear C-W equation in order to obtain the precision formation. For the formation flight examples cons idered above, the dynamics models are developed based on the leade r-follower (or hierarchical) structures. There are, however, some other ways to address formation system s. For example, Hadaegh and Goulet [65, 66] considered the thruster saturation a nd treated the formation system as a nonhierarchical structure. An adaptive satura ted controller is de signed based on this assumption. This adaptive controller con tinuously updates the estimated satellite mass and the mass is regarded as the uncertainty for this model.

PAGE 87

69 In this Chapter, the controllers based on the generalized model derived in Chapter 2 are designed. First, a performance criterion is discussed in Section 3.2. Then, in Section 3.3, in order to show the advantages of the sliding mode controller based on the full nonlinear model, controllers based on the following methodologies are develop: LQ method, synthesis, feedback linearization, and sliding mode controllers. Next, a comparison of the performance of these controllers is conducted in Section 3.4. The chapter concluded with a summary of the results obtained for the orbit control problem. 3.2 Performance Criterion Of The Leader-Follower Pair In this Section, the control performance criterion is given. In this Chapter, only the lower level leader-follower pairs control problem is addressed. Therefore, the criterion of the controller for the formation flying system can be simplified to the criterion of the controller for the leader-follower satellites pair. In most of the small satellite formation flying missions, such as UoSAT-12 and ACE, the relative position needs to be kept within a certain range and the attitude of the satellites needs to be tracked or maintained cooperatively. Therefore, the criterion of the controller is to minimize the above two branch errors (relative position and relative attitude errors) with minimum control efforts. Performance Criterion: The actual and desired relative position between the leader and follower satellites are ACTUAL and D ESIRED Then the error of the relative position is ERRORACTUAL DESIRED ]TAZ The error of the attitude is modeled as in the small Euler angle sense. Let the generated control, i.e., force [,ERROREL

PAGE 88

70 and torque commands, to be u Therefore, the criterion of the formation system is shown as Eq. (3.3), where 1Q 2Q and R are weighted matrices with proper dimensions. 0 D 121()2TTERRORERRORERRORERRORJQQu (3.3) TRudt In the following sections, linear and nonlinear controllers are developed based on minimization of the above criterion. Detailed performances for each method are explained in the corresponding section. 3.3 Relative Position Controller Design In this section, four control strategies are developed and implemented; they are LQR, feedback linearization, and sliding mode control. First, based on the linear dynamics model (C-W equation), the LQR and controllers are designed with uncertainties. Then, the feedback linearization method together with the synthesis is studied for the nonlinear model. Finally, the nonlinear control strategy, sliding mode control, is designed. All the controllers are assumed to have the full state feedback. The simulation results, discussion, and comparison are given in Section 3.4. 3.3.1 Linear-Quadratic State Feedback Regulator (LQR) The LQR method is used to design the optimal controller for the linear system in H 2 sense. 3.3.1.1 LQR Theory 3.1 ([67], p273): Assume the linear state space system is modeled as Eq. (3.4), where A B C and are constant matrices. The performance criterion is shown in Eq. (3.5) where 0Q and 0R are the weighted matrices with proper dimensions.

PAGE 89

71 x AxBuyCxDu (3.4) 01)2TTJxQxuRu dt (3.5) If matrix pairs (,)AB is controllable and (,)AC is observable, then there is one and only one optimal controller as shown in Eq.(3.6), in which P satisfies the Riccati equation (3.7). 1()TuRBPxt (3.6) 10TTPAAPPBRBPQ (3.7) 3.3.1.2 LQR synthesis model For the linear circular reference orbit, the relative position dynamics model can be written as the C-W equation as shown in Eq. (3.8), where a and 3Lnfa are the semi-major, gravitational parameter of the central body (Earth), and mean motion of the leaders orbit respectively. The state vector includes the relative position and velocity expressed in LVLH as [,,,,,]Tn x xyzxyz where subscript denotes the nominal model. The actuator, which is assumed to be n 11s is not easy to be synthesized in the LQR controllers design, therefore, the system equation is rewritten to include the actuator model before the design. First, the actuators transfer function is written in the system equation form as in Eq. (3.9), where the subscript represents the actuator. After synthesis of the actuator system equation with the C-W model as shown in Figure 3.1, the a

PAGE 90

72 combined linear system equation is as Eqs. (3.10) and (3.11), where the state vector is [,]Tna x xx 2200010000000001000000000100030002010000020001000000001n nn x xnnnn ua (3.8) aaaa x xuyx (3.9) Figure 3. 1 Combined state space 2200010000000000001000000000000100000030002010000000020001000000000001000000000100100000000010010000000001001nn x xunnAxBu (3.10)

PAGE 91

73 100000000000010000000000001000000000000100000000000010000000000001000000yxCxDu u (3.11) The controllers closed loop objective is to track and maintain the desired states. As such, 123456[,,,,,]Trrrrrrr is assumed to be the commands or desired states, and the error 61nerx is aimed to zero. In order to design the controller based on LQR synthesis, a new state is chosen as in Eq. (3.12), and new system equation is derived as Eq. (3.13). [,]T x xedt (3.12) 96966663660000AB x xuCI r (3.13) In order to minimize the errors that are last three elements of the state, weight functions for the state and control are chosen by try and error as in Eqs. (3.14) and (3.15). Using MATLAB LQR command, the feedback gain for the system can be obtained as in Eq. (3.16), where (:,:)ijK represents the column vectors from columns to i j NK AK and INTK are the gains for nominal model, actuator model, and integrator model respectively. 1212123312330.03500IQI (3.14)

PAGE 92

74 100010001R (3.15) (:,1:6)(:,7:9)(:,10:15)[,,][,,NAINTKKKKKKK ] (3.16) 3.3.1.3 LQR analysis model Figure 3.2 gives the analysis model of the LQR controller. The input is the desired relative position (i.e., commands), and the outputs are the relative position and velocity in LVLH coordinates and actuator forces. P and A V are the nominal C-W and actuator models. The integral model 1 s is used to reduce the steady state error. Furthermore, the uncertainty model is shown in Figure 3.3. P represents the C-W model with a uncertainty in the leader orbits mean motion that has 200% uncertainty. The un-modeled dynamics related to the nominal model and actuator are bounded by and where 50% and 20% uncertainties are assumed. 1 2

PAGE 93

75 Figure 3. 2 Analysis model of LQR controller without uncertainty Figure 3. 3 Analysis model of LQR controller with uncertainty

PAGE 94

76 Figure 3. 4 Typical loop gain and phase margin model In the classical control design, the gain and phase margins are used to study the robustness. The gain margin is the amount of gain increase required to make the loop gain unity at the frequency where the phase angle is 0 The phase margin is the difference between the phase of the response and 0 when the loop gain is 0 dB. It is generally found that gain margins of six or more combined with phase margins of 45 0 or more result in reasonable trade-offs between bandwidth and stability [46]. Each loop in the analysis model needs to be tested, and for simplicity, only one loop gain and phase margins model diagram is shown as in Figure 3.4. In this model, the loop gain and phase margins from command the to the output, relative position 1r x is tested. The control gains, time response and robust analysis based on LQR will be shown in Section 3.4. 3.3.2 /H Controller 3.3.2.1 Structured singular value ( ) theory /H is a tool used to design robust controllers and to analyze the robustness and performance of a system subjected to uncertainties, especially the structured

PAGE 95

77 uncertainties. The structured uncertainty represents the normalized uncertainty, which has a fixed structure with scale parts and matrix parts as shown in Eq. (3.17). i is a scale, ri I is a identity matrix, and irr i i is a square matrix. The /H controller searches only on the structured uncertainty, so it is not as conservative as based on the small gain theory. is defined as the inverse of the largest uncertainty, which causes the system to be unstable as shown in Eq. (3.18), where () is the singular value of the uncertainty and M is the open loop model. If the value of the closed loop system is less than 1, the controller is robust with respect to the structured uncertainties. 11[1, 12,,,j ]: ( 0} 22{[,,,,,,,11,[1,],],1}jjrrsrsFmmijidiagIIIisjF (3.17) 1min{):det()IM (3.18) synthesis is very difficult to be implemented, and D-K iteration is normally used [68, 69]. The basic ideas of the D-K iteration are first scaling the open loop plant by a state-space then compute a optimal D H control gain for the scaled plant. Second, compute set of optimal scaling for the closed-loop system under the fixed to reduce the K K bounds. Iterating these two steps until the peak value does not change much. The detailed information of synthesis theory and D-K iteration can be found in Lind [68] and Matlab [69]. 3.3.2.2 synthesis model First, the parametric uncertainty model, which has 200% uncertainty in the mean motion of the leader, is derived. Two correlated parametric uncertainties are considered.

PAGE 96

78 One is the mean motion of the leader cc 22022(1)W 11011(1)ccW where c. The other is the square value of the mean motion 20 n where cn 2220c 10 and 1 2 are normalized uncertainties. WW 24 12 ,i are uncertainties magnitude in this dissertation. Substituting these two parametric uncertainties into Eq. (3.8), the system equations can be rewritten as Eqs. (3.19) through (3.21). Therefore, the open loop parametric uncertainty state space model can be written as Eq. (3.22), which has the input-output relations as shown in Figure 3.5. [1,4]i are the exogenous signals, which includes disturbance and noise, ,[1,3iui ] are the control inputs, ,izi [1,4] are the outputs related to i and output ,[1,iyi 6] are composed of the relative position and velocity. 4101205110111202251012051112210120511232323232 x cxcxucWxcWxcxcxuZZcxcxu (3.19) 52042201204231204232222 14 x cxucWxcxuZcxu (3.20) 61033102210334210334 3 x cxucWxcxuZcxu (3.21)

PAGE 97

79 102020101012022011020001000000000000010000000000000100000030002011001000020000100100000000010013000000000000000020000000000020000000000000000000001000000000000010000000000000100000cccccWcWPcWcW 00 66671071060000000010000000000000100000000000001000000ABCD 0 (3.22) Figure 3. 5 Open loop model with parametric uncertainty The synthesis model is shown in Figure 3.6, where 1212{,,,}diag is the parametric uncertainty. M and A are normalized un-modeled uncertainties associated with the C-W equation and the actuator, respectively. M W and W are the A

PAGE 98

80 uncertainties magnitude with values of 50% and 20%. The measurement noise is assumed to be 660.02n W. The weighting function for the performance and control efforts are selected as I 660.10.1Perfss WI and 330.1(1)10uss W I respectively. The actuator is assumed to be the same as in LQR synthesis, which is 11s Figure 3.7 gives the input-output relations of this synthesis model. The differences between the synthesis model and the open-loop model are commands r noise n performance e, and control effort u. The feedback gain K between u and y are found through D-K iteration commands in MATLAB. Figure 3. 6 Synthesis Model

PAGE 99

81 Figure 3. 7 Synthesis model input-output relation 3.3.2.3 analysis model Figure 3. 8 Analysis model for nominal system Figure 3. 9 Analysis model for actual system After using the synthesis model to design the controller, the analysis model is constructed to evaluate the controller for the time response, frequency response, and robustness analysis, respectively. The analysis model for the nominal system and actual system are shown in Figure 3.8 and Figure 3.9. These two diagrams are similar to the ones in LQR analysis model and not discussed here. 3.3.2.4 NP, RS, and RP In /H type controller design, NP, RS, and RP are used to study the robustness and performance. RS: robust stability, family of models with uncertainty, which is stable.

PAGE 100

82 NP: nominal performance, family of models without uncertainty, which achieves the desired goals. RP: robust performance, family of models with uncertainty, which achieves desired goals [68]. After obtaining the feedback gain K, the closed-loop model with uncertainty is formed as shown in Figure 3.10. Furthermore, the equivalent models can be obtained as in Figure 3.11 and Figure 3.12. Here, 11122122(,)LSSSFSS GK is the lower Linear Fractional Transformations (LFT) of the models G and K with proper dimensions corresponding to the input and output. Figure 3. 10 Closed-loop model with uncertainty Figure 3. 11 Equivalent model 1 Figure 3. 12 Equivalent model 2

PAGE 101

83 In Figure 3.12, 11(S ) or 11S measures the robust stability, 22(S ) or 22S calculates the nominal performance, and ()S or S is the robust performance. If RS>1 and NP<1, the nominal system is easy to control, but the system may not be stable if the model changes. If NP>1 and RS<1, the system is easier to be stabilized, even if it varies by but a good performance is difficult to achieve. If NP<1 and RS<1, but RP<1, it is easier to stabilized for nominal system and also the nominal system is easier to control. But the good performance is difficult to achieve when the plant varies by One thing needs to be mentioned is the controller's order. Normally, the controller's order obtained through D-K iteration is high. If the order is larger than ten, it is difficult to implement for the real physical system. Therefore, Hankel model reduction is used here to reduce the controller's order. The detailed information of Hankel model reduction can be found in Reference [68] and [69]. The simulation results are discussed in Section 3.4. 3.3.3 Feedback Linearization The above two control strategies are mainly for the linear model. In the next, the generalized nonlinear dynamics model in Eq. (2.27) will be based and the partial or fully nonlinear control methods are used. 3.3.3.1 Feedback linearization [Slotine, 70, p207] The basic idea of the feedback linearization method is to algebraically map a nonlinear state space system dynamics into a (fully or partial) linear one so that the linear optimal control methodology can be applied. Feedback Linearization techniques can be viewed as ways of transforming original system models into equivalent models of a simpler form.

PAGE 102

84 Lie derivative is used in the feedback linearization and explained in appendix B. For the convenience, the dynamics system of the leader follower pair in the formation flight system is written as Eq. (3.23). [1,6](,)[(,)]LiLifxffxf is shown in Eq. (2.26) and (,)Lhxf is the output function as Eq.(3.24), where L f is the leaders true anomaly. Vectors ,1,2, 3 igi are listed in Eq. (3.25). '12122()(,)(,)LLL 2 x ffxfguguguyhxf (3.23) 123456(,)[,,,,,]TLhxfxxxxxx (3.24) 123[0,0,0,1/(1),0,0][0,0,0,0,1/(1),0][0,0,0,0,0,1/(1)]TTTgsgsgs (3.25) Rewrite Eq. (3.25) in the matrix form as 331415163324252634353600000000000000000001/(1)0001/(1)0001/(1)ggggGsgggsgggs (3.26) Constructing the mapping from x to () x and under this new state a fully/partial linear model appears. Now, the main steps of feedback linearization are explained. Step one is to calculate the relative degree for each output, where the relative degree is the minimum number of Lie differentiations of the output to have the control appear explicitly.

PAGE 103

85 Output1 : 111(,)()LLyhxfxf The first lie derivative of the output 1 is: 6'1111''14()()({1,0,0,0,0,0}fiiihyLhxxfxxxx )L Because the control does not explicitly appear, the second lie derivative is need. 21161'1''44141152()()()[()](){0,0,0,1,0,0}(,)ffiLiiLytLhxLhxxfxxx 163 f xfgugugu Therefore, the relative degree of the output 1 is r 12 (3.27) Output2 : 222(,)()LLyhxfxf 6'2221''25()()({0,1,0,0,0,0}fiiihyLhxxfxxxx )L 62'21''55241252[()]()(){0,0,0,0,1,0}(,)fiLiiLLhxytxfxxx 263 f xfgugugu

PAGE 104

86 Therefore, the relative degree of the output 2 is also r 22 (3.28) Following the same procedure for all the other four outputs, the relative degree vector is derived and shown in Eq. (3.29). 123456222111rrrrrrr (3.29) Step two is to change the original state space to a new one. A vector is defined according to the lie derivatives as Eq. (3.30). The governing equation for this new state space can be written as Eq. (3.31). 1,111,212,122,223,133,234,145,1566,1fffhLhhLhhLhhhh (3.30)

PAGE 105

87 1,21,11,22,22,12,23,23,13,2211,24222,2523,2634,14455,1566,16(,)(,)(,)fLfLLfffffffLhfxfLhfxffxfLhLhLhLhLhLhLh114151622425263343536000000000000000000uggguggguggg (3.31) Step three is to reduce the size of the set of states from step two. Seen in Eq. (3.31) 1,1 1,2 3,1 1,2 2,2 and 3,2 have explicitly show all the state variables in x therefore, the following simplified transform system is chosen as the new state space system. 1,11,22,12,213,13,221,2414151632,252425263,26343536000000000(,)(,)(,)LLLuufxfgggufxfgggfxfggg (3.32) Then the control force is designed to be Eq. (3.33). 411152336(,)()(,)()()(,)LLLLLLfxfuvfxfuGuvfxf 2fvff (3.33) Substituting Eq. (3.33) into Eq. (3.32), the linear model is derived as shown in Eq. (3.34). In this particular case, the feedback linearization gives the same results as the inverse dynamics, which is not the normal case.

PAGE 106

88 (3.34) 1,11,12,12,113,13,121,21,232,22,23,23,2000100000000010000000001000000000100000000010000000001vvv ] The linear model obtained above is not desirable since all the eigenvalues of the unforced internal dynamics are on the imaginary axis and the system is not asymptotically stable. In this case, most of the control efforts are wasted to stabilize the internal system. Therefore, a new feedback linearization control laws is proposed as shown in Eq. (3.35) with and Let e, the controller can be simplified as Eq. (3.36). Substituting this controller to Eq. (3.32), and defining a new state space as shown in Eq. (3.37), the final linearized state space system is obtained as (3.38). In the simulation, and are both chosen to be one. 0iaia 0ib ,1,1,[1,2,3diiii ib 1,24111,11,111,21,2112,25222,12,122,22,2233,13,133,23,23363,2(,)()()(,)()()()()(,)dddLdddLdddLfxfuabfxfuGababufxf vvv (3.35) 1,2411111112,2522222233333363,2(,)(,)(,)dLdLdLfxfuaebevfxfuGaeaebevufxf bev (3.36) 123123123456[,,,,,][,,,,,]Teeeeeeeeeeee T e (3.37)

PAGE 107

89 (3.38) 112213324435566000100000000010000000001000100100100010010010001001001eeeeveeveeveeee 3.3.3.2 Feedback linearization analysis model Based on the linear model (3.38) and the controller (3.36) derived from the feedback linearization, the controller architecture is shown in Figure 3. 13. The feedback gain K is obtained by the linear control theory, such as LQR or synthesis for the linear dynamics model. The nonlinear model is the relative position model as shown in Eq. (2.26). The closed-loop system for the linear model is shown in Figure 3. 14, where 123456[,,,,,]Teeeeee denotes the new state space. A is the same actuator model as before, and P is the linear model as shown in Eq.(3.38). The simulation results for the feedback linearization together with synthesis are shown in Section 3.4. Figure 3. 13 Closed-loop analysis model based on the feedback linearization

PAGE 108

90 Figure 3. 14 Closed-loop analysis model for the linear system 3.3.4 Sliding Mode Controller (SMC) The sliding mode control, a simple nonlinear control methodology based on the Lyapunov theory, is used in this section. A detailed discussion of the sliding mode methodology can be found in Ref. [70]. Sliding mode control Theory 3.2: Multiple-input Sliding Mode Control (Slotine, [70], p303) A nonlinear multi-input system is shown in Eq. (3.39), where represents the derivative with respect to time, )(in thin ijmmBb and 1[]imff Here f and B are denoted as state and input matrices. Also, B is assumed to be a square matrix, although some methods exist when it is not square in the multiple inputs sliding mode control. ()(1)(1)(1)(1)1(,,...,)(,,...,)iimnniijj inj x fxxxbxxx u mi,...,1 1,...,j m (3.39) In order to use the sliding mode control, two assumptions are needed to be satisfied. First, the matching conditions must be verified as shown in Eqs. (3.40) and (3.41), where i f and B are the predicted model of the state and input matrices. and iF

PAGE 109

91 ijD set the boundaries for the systems uncertainties. I is the identity matrix with a proper dimension. The second assumption is that the actual and predicted input matrices B and B must be invertible. If B and B are not invertible, the pseudo invertible is a common way in the sliding mode control. But in this dissertation, when B and B are not square or invertible, it is assumed that some directions thruster dont work and are detected as an error case i nr iii f fF 1,...,i m (3.40) () B IB ijijD i m,...,1 1,..., j m (3.41) Using subscript represents the desired states, based on these two assumptions, the sliding surface is defined to be Eq. (3.42), where d d x xx is the tracking error and must be strictly positive. The optimal control is shown in Eq. (3.44), where and 12(1)(1)(1)(1)12mTnnnrrrmxxxx (1inrix ) is defined by Eq. (3.43) and Controllers parameter should be chosen according to Eq. (3.45), where 1sgn()[sgn()]iiksks0i m ik based on the stability condition (1)(1)(1)()iinniiiidsxxdt inrix (3.42) (1)(1)1iiinnnridiii x xx (3.43) 1(1)[sgn()]nruBxfks (3.44) (1)1)inniiiijjiijriijjijDkDkFDxf (1 (3.45)

PAGE 110

92 Sliding mode control has strong relations with the classical PID controller, therefore, this controller is easier to understand and convenient for online tuning in implementation. If we rewrite out system model to be Eq. (3.46), and substituting Eq. (3.44) into it, Eq. (3.47) is arrived. Using Eq. (3.43), the final relation is shown in Eq. (3.48). ()(1)(2)(1)(,,,...,)n n x fxxxxBu (3.46) ()(1)(2)(1)1(1)111(1)111(1)(,,,...,)[sgn()]sgn()sgn()nnnrnrnrxfxxxxBufBufBBxfks f BBfBBxBBksfBBfBBksBBx (3.47) 111111()(1)(1)sgn()nnnnd x fBBfBBksBBxBBxBBxdt (3.48) In this expression, 1(1)n B B x represents the proportional gain, 1(1)n B Bx and 11nd B Bdt x are similar to the differential gain, 1 f BBf accounts the uncertainties, and 1sgn() B Bks relates to the thruster bandwidth and how faster the thruster can force the performance switch between the sliding surfaces. 3.4 Relative Position Simulation Results and Discussion Initially, the leader is in a 500 Km altitude circular equatorial orbit with a true longitude at epoch () of The follower is also in a 500 Km circular orbit, but has a inclination and a true longitude at epoch (u) to be The small inclination and true longitude at epoch difference are chosen to keep the relative distance between 0Lu 00 00.1 0F 00.1

PAGE 111

93 the leader and follower satellites on the order of 10 Km. The thruster is assumed to work at 10 Hz in this section. The four controllers designed in the previous sections are tested and the simulation results are compared. 3.4.1 State Feedback Gain Here, the feedback gain for LQR, feedback linearization, and sliding mode controllers are listed. LQR Controller (Based on C-W equation): 1.4160.00202.5930.00100.0021.41600.0012.5930001.416002.593NK (3.49) 1.4940001.4940001.49AK 4 (3.50) 0.1870000000.1870000000.187000INTK (3.51) Controller(Based on C-W equation): MMMMMABPCD (3.52) 32.2750.0020.0010.0060.0010.0010.0020.0790.02900.03200.00100.0030.01100.001000.0220.0090.0020.00200.0010.0010000.0020.0020.0010.001001100000.0070.0060.0010.0010.0010000000.0060.0MA0700000000.0020.0020000000000.0020.00600000000.0030.006 (3.53)

PAGE 112

94 61.58410.8545.43391.91115.4159.9681.3033.8114.4041.7066.9858.3030.4183.2893.3020.3506.4125.4630.1361.5221.1360.0542.8391.9040.2121.5361.4330.2133.0332.6320.2770.8391.5950.3971.582MB3.0480.3031.4841.3990.3282.6152.5250.4850.3540.0510.4860.5670.0952.6510.2690.1373.9770.3300.242 (3.54) 21.1070.0210.0060.0010.0020.0050.0050.0020.046100.1850.0790.0710.0350.0330.0170.0310.0060.0040.1100.0940.0620.0250.0290.0340.0300.0020.003M C (3.55) 0.0110.0010.0010.0260.0060.0020.001000.002000.001000.00100MD (3.56) Controller in Feedback Linearization (Based on the linear model obtained by feedback linearization): FMFMFMFMFMABPCD (3.57) 2.061002.0260.0050.05902.06100.0062.0270.026002.0610.0590.0262.0260000.5040000000.5040000000.5FMA 04 (3.58) 30.2880.2760.6120.0330.0320.0700.5850.2230.3760.0670.0260.0430.3280.6380.1330.0380.0730.015100.1570.1650.3390.0130.0140.0280.3290.1190.2100.0270.0100.0170.1840.3540.0870.0150.0290.007MB (3.59) 30.3040.6180.3470.1270.2680.150100.2910.2350.6740.1340.0970.2880.6460.3970.1410.2760.1710.071M C (3.60)

PAGE 113

95 80.046000.161001000.046000.1610000.046000.161MD (3.61) Parameters in SMC (Based on the generalized nonlinear model): [100,100,520]TK (3.62) [0.1,0.1,0.1]T (3.63) 3.4.2 Robustness and Performance Analysis The robustness and performance are compared for the first three controllers, which are LQR, and feedback linearization plus Table 31 Gain and phase margins Table 31 Gain Margin (dB) Phase Margin 0 LQR 6.15 77.06 8.19 77.16 gives the gain and phase margins for the LQR and synthesis, for both methods. The controllers designed are reasonable in the robust, performance, and stability sense for the gain margin is larger than 6 dB, and the phase margin is greater than 45 0 lists the robustness and performance results for the and feedback linearization controller. NP, RS, and RP are all less than one, which means that the system can achieve the desired performance even if the system has the predefined structured uncertainties. The controllers order is 18 for both cases and difficult for implementation in real physical system, and in the simulation, Hankel model reduction method is used. After the model reduction, controllers orders are 9 and 6 separately for and feedback linearization. Table 32

PAGE 114

96 Table 32 Robustness performance and controller order comparison NP RS RP Controller Order before Reduction Controller Order after Reduction 0.57 0.69 0.84 18 9 Feedback Linearization Plus 0.2 0.734 0.796 18 6 3.4.3 Time Response Based on C-W Model Figure 3.15 through Figure 3.20 show the relative position maintaining performance of the above four controllers for the C-W equation. Due to the much longer settling time of the feedback linearization plus controller, its simulation results are shown individually. From Figure 3.15, Figure 3.17, and Figure 3.19, it can been seen that LQR controller has much longer rising time and larger overshot than and SMC. Furthermore, the controller gives better performance than SMC in rising time and overshot senses. Therefore, in the time response, and SMC controllers are the good choices for this simulation. Figure 3. 15 Relative position error in x Figure 3. 16 Relative position error in x

PAGE 115

97 Figure 3. 17 Relative position error in y Figure 3. 18 Relative position error in y Figure 3. 19 Relative position error in z Figure 3. 20 Relative position error in z Figure 3.21 through Figure 3.22 give the control effort in and SMC controllers. In this sections simulation, the thruster has no maximum limitation. The stable state of the controller in x direction is not zero due to the fixed steady state error in x direction, which is not shown clearly in Figure 3.15. In Table 33, the fuel consumption is calculated by the norm 2 value of the control action for the first three hundred seconds. SMC saves up to 30% as that of controller. The SMC gives lots of chattering, which will be discussed later in this Chapter.

PAGE 116

98 Figure 3. 21 Control force in x Figure 3. 22 Control force in y Figure 3. 23 Control force in z Table 33 Fuel consumption Fuel Cost (KN/Kg) 0.0492 SMC 0.0362 3.4.4 Time Response Based on C-W Model with Uncertainties In this section, the uncertainties in the parametric and un-modeled dynamics are considered. The mean motion of the leader is assumed to 200% of the exact value. The un-modeled dynamics is chosen to be 1.5 times of the exact value, and the actuators magnitude is 20% above the original one. Same reason as the simulation with the exact model, the feedback linearization method gives the much longer settling time, and the results are shown separately in Figure 3.25, Figure 3.27, and Figure 3.29. In y direction, a

PAGE 117

99 fixed stable state error exists, with magnitudes to be 2 Km. The LQR, and SMC give approximately the same stable state error, whereas, the LQR has the largest overshot, and longest settling time in x and y direction. These simulation results demonstrate that the control is stable under the maximum uncertainties condition, and the performance is degraded by a fixed error. When the simulation is zoomed in, the chattering phenomena of the SMC are still there. Figure 3. 24 Relative position error in x Figure 3. 25 Relative position error in x Figure 3. 26 Relative position error in y Figure 3. 27 Relative position error in y

PAGE 118

100 Figure 3. 28 Relative position error in z Figure 3. 29 Relative position error in z The control efforts of the and SMC are shown in Figure 3.30 through Figure 3.32. Under the uncertainties condition, the needs much larger efforts than SMC, and SMC saves up to 99% in the fuel consumption of the first 3000 seconds simulation. Figure 3. 30 Control force in x Figure 3. 31 Control force in y

PAGE 119

101 Figure 3. 32 Control force in z Table 34 Fuel consumption Fuel Cost (KN/Kg) 18.8846 SMC 0.1138 3.4.5 Time Response Based on Actual Nonlinear Model with Thruster Limitation The above two simulations have not considered the maximum force the thruster can provide. Also, not like feedback linearization and SMC, the linear controllers have not been implemented in the real nonlinear model. Therefore, in this section, the actuator is chosen to have the same properties as IS120 (ST-5). The maximum thruster force is 2 N. Because the parametric and un-modeled uncertainties in the dynamics and actuator models have been studied in Section 3.4.4, they are not included here. Based on the above two simulations, and SMC are chosen to be compared in this section. Figure 3. 33 Relative position error in x Figure 3. 34 Relative position error in y

PAGE 120

102 Figure 3. 35 Relative position error in z Figure 3.33 through Figure 3.35 show the relative position results. The controller is unstable, whereas the SMC has a stable results. The SMC is based on the nonlinear model and robust to the model uncertainties. However, is based on the linear model and robust to the bounded structured uncertainty. According to Chapter 2, the open loop relative position difference between the C-W equation and nonlinear generalized model is unbounded. Therefore, the controller is not able to give a stable results in the real nonlinear system under the chosen actuator. The fuel consumption of the SMC is 0.1143 KN/Kg of the norm 2 sense. 3.4.6 Conclusion Based on the four control methodologies designed and the simulation results shown above, the following conclusions can be obtained. The LQR, and feedback linearization plus controllers designs are a standard procedures with uncertainties considered formally. Nevertheless, for SMC, more parameters need to be tuned and the processing time may be longer and the uncertainties have no fixed form.

PAGE 121

103 The LQR, feedback linearization and SMC is easier to understand and calculated than synthesis. Also, SMC is similar to PID controller, which makes the online tuning possible. The LQR and synthesis is mainly used for the linear model, whereas, the feedback linearization and SMC are two of the main techniques for the nonlinear system and robust with respect to the unstructured nonlinearity. Furthermore, for the control problems in this dissertation, tuning the weights of the performance and control is very difficult. In SMC, the chattering phenomena exist. As shown in the simulation, the SMC costs less fuel than controller, which makes it helpful in the space mission. Therefore, Sliding mode control is chosen as the basic control method in the leader follower pairs maintaining and tracking problem. Up to now, the orbit part of the low level controller has been designed. Been compared with the LQR method, synthesis, and feedback linearization, the sliding mode controller is chosen as the controller in this dissertation.

PAGE 122

CHAPTER 4 CONTROL STRATEGIES: ATTITUDE AND ORBIT 4.1 Introduction In Chapter 3, after compared with LQR, /H synthesis, and feedback linearization, sliding mode control is chosen as the basic control methodology applied for the formation flying system in this dissertation. In this Chapter, first, the sliding mode controller is designed for attitude in Section 4.2. After that, the coupled relative position and attitude controller is developed in Section 4.3. Although the classical sliding mode controller performs reasonably well, it does result in actuator chatters. In Section 4.4, boundary layer control is proposed to mitigate the chattering phenomenon. A genetic algorithm is used to find the local sub-optimal solution in Section 4.5, finally. 4.2 Relative Attitude Control 4.2.1 Control Architecture According to the formation mission described in Chapter 2, with the assumption that both the leader and follower satellites are gravity gradient stabilized, an attitude controller is required to make the follower satellite point to the leaders nadir point. Therefore, the attitude controller is composed of an inner loop and an outer loop as shown in Figure 4.1. The inner loop controller tracks the relative desired attitude of collectively pointing assuming gravity gradient stability, whereas the outer loop controller maintains the gravity gradient pointing, which will not be considered in this dissertation. 104

PAGE 123

105 Figure 4.1 Attitude control architecture The absolute angular velocity of satellites are composed of o and B where o is the angular velocity for gravity gradient pointing, and B is the relative angular velocity associated with tracking the relative desired attitude. As shown in Chapter 2, satellites relative rotation sequence with respect to gravity gradient pointing is 1-2. Also, expressed in LVLH frame, ()0LVLH has only one nonzero quantity which is in direction. The relative angular moment of satellites is denoted by Eq.(4.1), then the time derivative of Eq. (4.1) is calculated by Eq. (4.2). Furthermore, based on the assumption that the inertia matrix of the satellites is diagonal and the relative desired attitude ( L f y 123,,T ) is small, the last term in Eq. (4.2) only has one nonzero element, which is in the third axis as shown in Eq. (4.3) and not concerned for the relative rotation from the gravitational gradient pointing attitude to the coupled desired attitude. Therefore, the relative attitude model is shown as in Eq. (4.4) and used in all the studies later. ()()BrelBHJ B (4.1)

PAGE 124

106 ()()()()()()()()()()()()(BBBBrelrelBoBBBBBoBBBBBBBoHJHJJJ )BBJJ (4.2) 321()()312213221111112231101110000BBoBLLLLLLJJfJJfJfJfJfJf 12J (4.3) ()()()()BBBBBBBJJ T (4.4) 4.2.2 The Inner Loop Controller In this section, a sliding mode controller is designed for the inner loop relative attitude maneuver. The designed sliding mode controller is demonstrated to be able to track and maintain the desired states under uncertainties. The attitude governing equation expressed in Euler angles form is shown in Eq. (2.77). The error angles, i.e., elevator and azimuth error angles are derived in Eqs. (2.69) and (2.72). The corresponding sliding mode controller is designed in the same way as shown for the relative position model in Section 3.3.4 Before the sliding mode controller is demonstrated to be capable of dealing with uncertainties. The first simulations are used to show the effects of uncertainties on the relative position and attitude model. If the relative position model is written as () x fxBu then modeling uncertainties considered here are established as

PAGE 125

107 (1)f f f and (1)b B B where f and b are state and input matrices magnitude uncertainties. Also the uncertainties related to the inertia matrix for each satellite are modeled as (1)inertiaJJ in which i nertia is the magnitude uncertainty for the inertia matrix. As listed in Table 4-1, these uncertainties are considered in different cases and compared through simulations. The difference between the exact model and model with uncertainties is calculated by iExacx t iUncerx where i x denotes the state. Just as used in Chapter 2, the inertia matrix of the satellite is assumed to be thi 200429 4430012900J Kgm The leader satellite is in a 500 Km altitude circular equatorial orbit with an initial true longitude (u) of 0. The follower satellite is also in a 500 Km circular orbit, but has an inclination and a true longitude () of The small inclination and true longitude were chosen to keep the relative distance between the leader and follower satellites on the order of 10 Km. oL o 0.1o oFu 0.1o 5% n ertial 5% n ertial fb Table 41 Cases Case 1 Exact Model Case 2 i Case 3 Without control i and 10% As shown in Figures 4.2 and 4.3, the differences between cases 1 and 3, are compared in one orbit period, the z direction difference goes up to 4 Km, which is approximatly 25% of the desired relative position. Furthermore, the difference is unstable and increases with time.

PAGE 126

108 Figure 4.2 Relative position difference between exact and uncertainties model 0 2000 4000 6000 -4 -2 0 2 4 6 Comparison: Exact & Uncertainties Model R e l at i ve pos i t i on difference in z directionsTime Seconds Figure 4.3 Relative position difference between exact and uncertainties model Under the same initial condition in angular velocity, Figures 4.4 and 4.5 show the inertia matrix uncertainty with the same magnitude of the satellite does not affect the attitude (solid line, case 1 and case 2) in the open-loop case because the inertia matrix chosen here is a diagonal matrix. In the attitude model 0JJ the same result is obtained if the inertia matrix J is the diagonal and the uncertainty is also in the diagonal entries. However, the uncertainties of the relative position model have effects on the attitude performance due to the coupling between orbit and attitude. The dotted lines in Figures 4.4 and 4.5 show the effects on the attitude due to the relative position models uncertainties. The elevation error angle is related to the z direction of the relative position. Seen from Figure 4.2 and Figure 4.3, in the range of 0Lf when the relative position in z direction of the uncertainty model is higher than that of the exact model, i.e., 0zExactzUncer the absolute value of the elevator error angle for the exact model is smaller than that of the uncertainty model. The rotation from the actual to the desired attitude position is first through an elevation error angle against the negative body

PAGE 127

109 axis 1, therefore, just as shown in Figure 4.4, case 1 minus case 3 is large than zero. The same reasoning can be used to explain Figure 4.5. Figure 4.4 Elevator error angle difference between the exact and uncertainties models Figure 4.5 Azimuth error angle difference between the exact and uncertainties models Figure 4.6 Elevator error angle if no control is applied Figure 4.7 Azimuth error angle if no control is applied Figure 4.6 and Figure 4.7 give the simulation results of the error angles for the exact mode without control efforts. In order to satisfy the attitude criterion, a controller needs to be designed to drive these two error angles to zero. In Table 4-2, all five cases use the same designed sliding mode control and have the same controller parameters.

PAGE 128

110 Different uncertainty cases are simulated. In the simulation, U denotes the uncertainties in fb f and and b J U represents uncertainty in inertial Table 4-2 Cases Case 1 Exact Model Case 2 5%inertial Case 3 10%fb Case 4 10%fb and 5%inertial Case 5 [50;80;0.01]TaK[5;2;0.01]Ta [5%,5%]inertial and ,[10%,10%]fb random Figure 4.8 and Figure 4.10 show that the control precisions for both elevation and azimuth error angles are almost the same in the fixed uncertainty case, which are on the order of O radians or arc seconds. But the effects of the different uncertainties can be seen from Figure 4.9. If the model uncertainties only exist in the relative position model, the results are almost the same as the exact model. However, if satellites inertia matrix has uncertainty, the settling time and overshot are larger than the case when the inertia matrix is exact. In these simulations, for example, the settling time for cases 1 and 3 is nearly 400 seconds, whereas that of case 2 and 4 is 500 seconds. This is because the controller designed based on the uncertainties inertia matrix is not exact for the model. Figure 4.11-Figure 4.12 show the results if uncertainties are random (Case 5). The sliding mode controller works well in this case. The peak values in Figure 4.11 happen when large uncertainties appear. 4(10) (10)O

PAGE 129

111 Figure 4.8 Elevator error angle under the sliding mode controller Figure 4.9 Azimuth error angle under the sliding mode controller Figure 4.10 Azimuth error angle under the sliding mode controller Figure 4.11 Elevator error angle in case 5

PAGE 130

112 Figure 4.12 Azimuth error angle in case 5 4.3 Combined Relative Position and Attitude Controller In the previous sections, controllers have been designed for the relative position and attitude model separately. Here, both the position and attitude controllers are turned on. Performance is studied to show whether or not the combined controller is better than the case when only one controller is working. Different simulated cases are listed in Table 4-3. The sliding model controllers parameters chosen for this section are , [50;80;0.01]TaK [5;2;0.01]Ta [100;100;520]TpK and Furthermore, as shown in Table 4-3, the simulation for cases 1 and 3 is running only for 1 period, whereas 22 periods is simulated for case 2 in order to get a good understanding for a longer simulation. Then based on case 3, the fuel consumptions for the force and torque commands are calculated. In this section, the thruster provides a 2N maximum force command (the same characteristic as IS120) and the maximum torque command that can be supported is 12 mN-m based on TW-4A12. [0.1;0.1;0.1]Tp

PAGE 131

113 Table 4-3 Cases Case 1 1 period simulation Exact Model Case 2 22 periods simulation Exact Model Case 3 [50;80;0.01]TaK[5;2;0.01]Ta [100;100;520]TpK [0.1;0.1;0.1]Tp 1 period simulation [5%,5%]inertial ,[10%,10%fb and ] random Figure 4.13 Figure 4.15 give the relative position results for the simulation cases 1 and 3. Because the attitude model does not affect the relative position model, whether the attitude controller is on or off does not influence the relative position performance. The simulation shows that the sliding mode controller rejects the uncertainty. Figure 4.16 Figure 4.19 show the performance of EL and AZ, i.e., elevation and azimuth error angles, when both controllers are on. The precision, which is on the order of arc seconds, is almost the same as the case when only the attitude control is on. Therefore, from this point of view, the combined controller is not very necessary. The purpose of the simulations conducted here is to show that the sliding mode controller has the capability of maintaining both the attitude and position. (10) Figure 4.13 Relative position in x direction Figure 4.14 Relative position in y direction

PAGE 132

114 Figure 4.15 Relative position in z direction Figure 4.16 The elevation error angle Figure 4.17 Zoomed in elevation error angle Figure 4.18 Zoomed in azimuth error angle Figure 4.19 Zoomed in azimuth error angle Figure 4.20 Zoomed in force commands in x direction

PAGE 133

115 Figure 4.21 Zoomed in force commands in x direction Figure 4.22 Force commands in y direction Figure 4.23 Zoomed in force commands in z direction Figure 4.24 Torque commands in axis 1

PAGE 134

116 Figure 4.25 Zoomed in torque commands in axis 2 Since the state performances of case 2 are similar to case 1, only the force and torque commands results are shown here in Figure 4.20 Figure 4.25. The fuel consumption is calculated by 2-norm of the force or torque commands. In this case, the total fuel consumption for the force commands is 0.043, and that of the torque is 17.815 With the exception the x direction force command, all the other control commands switch between their maximum and minimum values. This type of actuator switching is energy inefficiency, so methods are pursued to reduce this waste. The approach is to appropriately select the sliding mode control parameters and that is the subject of the remainder of this chapter. /KNKg NM 4.4 Boundary Layer Control In all the previous sections, the sliding mode controller is shown the capability of tracking and maintaining the desired states. However, there are two main disadvantages for the sliding mode controller. One is the chattering phenomenon, which is explained in Section 4.4.1 and the other is that performance has not been optimized, which will be

PAGE 135

117 discussed in Section 4.5. In this section, the boundary layer cont rol or soft switching control is introduced to reduce the chatters. 4.4.1 Chattering Phenomenon In previous sections, which discussed the sliding mode control, several disadvantages are demonstrated, such as lo w precision, chattering phenomenon and high fuel consumption. As stated in Ref. [74] chattering phenomenon is generally perceived as a motion, which oscillates about the sliding manifold. Ref. [74] further stipulates that there are two reasons causing this phenomen on. One is the switching delay in real hardware, and the other is that the hard sw itching, for example, the signum function in sliding mode control, can cause the state osci llating in a high fre quency. Switching delay in realities may occur because of the hardware design or the parasitic dynamics in series with the plant. Chattering is problematic becaus e it will result in actuator failure and may also excite the unstructured uncertainties in the system ([74]). Furthermore, chattering might affect the stability and performance. 4.4.2 Boundary Layer Control Normally, the boundary layer controller is used to mitigate the chattering phenomenon. Boundary layer control uses th e so-called piecewise linear or smooth approximation of the switching element in the bo undary layer of the sliding manifold [75, 76, 77, 78]. The basic ideas can be shown in Figure 4.26 and Figure 4.27. Here, SMC denotes the sliding mode control and SMC + BLC represents the sliding mode plus boundary layer control. When the SMC method is used, the actual motion has many oscillations because of the complete relaxati on that exists in the region. But in SMC + BLC, the sliding surface includes a hard sw itching outside of the boundary, and a soft

PAGE 136

118 switching inside of the boundary. The performance of the BLC is smoother and less fuel consumption. Figure 4.26 Sliding surface in SMC Figure 4.27 Sliding surface in SMC + BLC In this section, the following boundary layer saturation function shown by Eq. (4.5) [79] is used to replace the signum function in the sliding mode control as shown in Eq. (3.44). iiiiiissssats)sgn()()sgn( (4.5) 11iiss 1,2,...,6i 4.4.3 Simulation Results For the same initial orbit information as used in Section 4.2.2, the following parameters are chosen for the sliding mode and boundary layer controller. The simulation results are showing detail for these relative position and attitude performance, sliding surface, and control efforts. Case 1 [100;100;520]TpK [0.1;0.1;0.1]Tp [50;20;0.01]TaK , [5;3;0.01]Tap [0 .01;0.01;0.01]T [0 .02;0.007;0.02]T a

PAGE 137

119 Figure 4.28 Relative position in x direction, transient stage. Figure 4.29 Relative position error in x direction, steady stage. Figure 4.30 Relative position in y direction, transient stage. Figure 4.31 Relative position error in y direction, steady stage. Figure 4.32 Relative position error in z direction, transient stage. Figure 4.33 Relative position error in z direction, steady stage.

PAGE 138

120 Figure 4.34 Sliding surface for the relative position in x direction, transient stage. Figure 4.35 Sliding surface for the relative position in x direction, steady stage. Figure 4.36 Sliding surface for the relative position in y direction, transient stage. Figure 4.37 Sliding surface for the relative position in y direction, steady stage. Figure 4.38 Sliding surface for the relative position in z direction, transient stage. Figure 4.39 Sliding surface for the relative position in z direction, steady stage.

PAGE 139

121 Figure 4.28 Figure 4.33 give the relative position. Compared with SMCs simulation, SMC + BLC mitigates the chattering phenomenon in the stable stage, and the error goes to zero which gives a high resolution. From Figure 4.34 Figure 4.39, sliding surface smoothly arrives the target surface, which is zero. Although sliding surfaces in x and y directions have oscillations, they are very small and on the order of and respectively. Seen in Figure 4.38, the error is outside of the switching region in the beginning, and the thruster commands go into saturation. In about 500 seconds, the error is within the desired region, then the boundary layer controller works to drive the error to zero along a smooth trajectory. 14(10) 13(10) Figure 4.40 Figure 4.43 show the transient and steady stages for the error angles, i.e., the elevation and azimuth error angles. In the sliding mode control only case, the precision is on the order of (1) arc seconds with chatters. Seen from the current figures (Figures 4.40-4.43), when boundary layer controller is working together with the sliding mode controller, the elevation error angles precision is on the order of radians or arc seconds. In the same time period (i.e., steady state), the azimuth error angle is zero and without oscillation. The sliding surfaces for the attitude are not shown here for they are similar to the sliding surfaces for the relative position. 11(10) 5(10)

PAGE 140

122 Figure 4.40 The elevator error angle in the transient stage Figure 4.41 The elevator error angle in the steady stage Figure 4.42 The azimuth error angle in the transient stage Figure 4.43 The azimuth error angle in the steady stage The fuel consumption and maximum force comparisons are also calculated. It is observed that SMC+BLC reduces the total fuel consumption by 79% of that of the SMC in the relative position controller. Similarity, for the attitude part, the fuel consumption is reduced by 90% of the case where only SMC is used. Furthermore, for the sliding mode plus boundary layer controller, the maximum force command required is only half of the case when SMC is applied alone.

PAGE 141

123 The force and torque commands are shown in Figure 4.44 Figure 4.53. The shape of these force commands is similar to the sliding surface for the relative position shown in Figure 4.34 Figure 4.39. Nevertheless, the disadvantage against the SMC is that the SMC + BLC has a larger settling time, which is nearly one orbit period in the attitude performance. Figure 4.44 Force commands in x direction of the local frame, in transient stage Figure 4.45 Force commands in x direction of the local frame, in steady stage Figure 4.46 Force commands in y direction of the local frame, in transient stage Figure 4.47 Force commands in y direction of the local frame, in steady stage

PAGE 142

124 Figure 4.48 Force commands in z direction of the local frame, in transient stage Figure 4.49 Force commands in z direction of the local frame, in steady stage Figure 4.50 Torque commands in body axis 1, in transient stage Figure 4.51 Torque commands in body axis 1, in steady stage Figure 4.52 Torque commands in body axis 2, in transient stage Figure 4.53 Torque commands in body axis 2, in steady stage

PAGE 143

125 In general, it has been demonstrated that while using the sliding mode controller together with the boundary layer control, the relative position and attitude can be maintained with small control commands. Nevertheless, one important point has not been addressed, which is how to choose the controllers parameters in order to obtain an optimal or a sub-optimal solution. In next section, a genetic algorithm will be proposed as an effective way to select the controllers parameters. 4.5 Genetic Algorithm In Section 4.4, the sliding mode together with a boundary layer controller was shown to provide a good performance (i.e., a high precision and less chattering phenomenon). In this section, we discuss how to choose appropriate parameters for this type of controllers. First, the optimization problem is stated, then a solution based on genetic algorithm is proposed. Finally, simulation results demonstrate that the genetic algorithm is able to find the local minimum or sub-minimum solutions. 4.5.1 Optimization Problem The system model is listed as in Eq. (4.7), where the state vector is composed of the relative position and Euler error angles as shown in Eq. (4.6). The sliding mode control (SMC) and the boundary layer control (BLC) are used together to maintain the system in a real time desired state as in Eq. (4.8). The sliding mode controller is designed to be Eq. (4.9), where the predicted models B and f must satisfy the conditions whosn in Eqs. (4.10) and (4.11). 61123[,,,,,]Txyzx (4.6) (,,) x fxxtBu (4.7)

PAGE 144

126 ()DESIRED x gx (4.8) 1[rsuBxfKsat ()] (4.9) ,1,2,...,iii f fFi n (4.10) () B IB ;,1,2,...,ijijDijn (4.11) In Eq. (4.9), riDESIREDiii x x x and the state error x is calculated by D ESIREDxx For stability, 0i is required. Also, (s )Ksat is defined as 6616...,()]TsKsat 11(),sKsat ()[sKsat where 61:[]PositionAttitudeiKKKK ; and 0:i 61;]tionAttitude :[iPosi ...,n The restriction condition on parameter is shown in Eq. (4.12). The sliding surface is defined to be Eq. (4.13). The saturation function in the boundary layer is shown as Eq. (4.14). Therefore, in this optimization problem, total 18 parameters need to be optimized as listed in Eqs. (4.15) and (4.16). The performance criterion is in a quadratic form as shown in Eq. (4.17). 1,2,iKi 1(1)niiiijjiijrijjijDkDkFDxf i where 0i (4.12) ~()iiiidsxxdt rix (4.13) 1||1()1,2,...6||1iiiiiississ sat (4.14) In SMC: 31313131PositionPositionAttitudeAttitudeKK (4.15)

PAGE 145

127 In BLC: 3131PositionAttitude (4.16) Performance Criterion: 01[()()]2tTTDESIREDDESIREDJxxQxxuRu dt (4.17) 4.5.2 Genetic Algorithm In the sliding mode and boundary layer control method, there are several disadvantages. For example, it may take a long time to find a set of appropriate parameters for the controller. Furthermore, there is no guarantee that the parameter set selected is the optimal or sub-optimal solution [80, 81]. In designing the sliding mode and boundary layer control, a total of 18 parameters are required. Also, the system is not monotonic and the evaluation takes a long time. Therefore, the conventional optimal method, such as the large-scale optimization given by Powell [82], does not work efficiently. In this dissertation, the genetic algorithm is used to find the local minimum solution for the sliding mode and boundary layer controller automatically. The genetic algorithm performs a parallel, non-comprehensive search for the global maximum/minimum of the graph [83]. The search is weak because there is no guarantee that the global maximum/minimum will be found. However, the result should be a good approximation value. If there exists a good specialized optimization method for a specific problem, then genetic algorithm may not be the best optimization tool for that application. Most users of genetic algorithms typically are concerned with problems that are nonlinear. This often implies that it is not possible to treat each parameter as an independent variable, which can be solved in isolation from the other variables. There are six parameters in the genetic algorithm; these are the number of parameters need to be optimized, the bit length, the population size, the number of

PAGE 146

128 generation, the crossover and mutation probability. Bit length is how many bits are used to represent a chromosome or parameter. The real value for each parameter or chromosome bl x is calculated by 11minmaxmin2()21blkkkblg x xxmin x where is the i gene in the chromosome [82, 83]. ig th x and max x define the range for the parameter x The population size represents the number of chromosomes in one generation. The number of generation determines how many generations the simulation has. The crossover probability describes the chance that bits of the same position in different chromosomes exchange their values, whereas, the mutation probability denote the chance bits in different position of the same chromosome exchange their values. The process of evaluation, selection, crossover and mutation forms one generation in the execution of a genetic algorithm as shown in Figure 4.54 through Figure 4.57. Figure 4.54 A generation of the chromosome First, evaluation of the objective function based on the current generation (as shown in Figure 4.54) is calculated. The fitness function is obtained from these evaluation values. The fitness function transforms the measurement of performance into an allocation of reproductive opportunities [82]. The evaluation of a chromosome i f is independent of the evaluation of any other chromosomes. However, the fitness of that

PAGE 147

129 chromosome is always defined with respect to other members of the current population. Typically, fitness is defined by /ii f f where i f is the average evaluation of all the chromosomes in the current generation as Figure 4.55. Figure 4.55 Evaluation, fitness, and selection Figure 4.55 also shows the selection step (the triangle shapes above the body), which determine whether a particular chromosomes gene will stay in the next generation. In canonical genetic algorithm the remainder stochastic sampling method is used based on the fitness of each chromosome. In this method, the integer part of the fitness is the number of copies will be used in the next generation, whereas the remainder part is compared with other chromosomes remainder part to determine who goes into the next generation. For example, if the population size is 4 and the current fitness values are 1.3, 1.1, 0.6 and 1.7. Chromosome 1, 2, and 4 must have one copy in the next step because the integer part is 1. The remainders, 0.3, 0.1, 0.6 and 0.7, determine that the chromosome 4 has another copy inherited by the next generation. For remainder part of the chromosome 4, 0.7, is larger than all that of the other chromosomes, the chromosome 4 will have two copies in the next step.

PAGE 148

130 Figure 4.56 Crossover After selection has been carried out, the crossover is applied to randomly paired chromosomes with a probability denotes as shown in Figure 4.56. For each bit in the chromosome, (i) determines whether this bit is going to be exchanged with the corresponding bit in the other chromosome. This is so called uniform type crossover. Typically, has a high probability. cP ciP 1,2,...,ciP bl Figure 4.57 Mutation Seen in Figure 4.57, mutation is the last step before the next generation is generated. Probability of the mutation is normally a low value. The operation will randomly exchange the bits within the chromosome.

PAGE 149

131 After the process of selection, crossover, and mutation is complete, the next population can be evaluated. The basic idea of this procedure is that the good genes will have higher chance to be inherited by the next generation. 4.5.3 Simulation Results For the results shown below, the optimization was done over one orbital period. The results from the previous SMC+BLC without genetic algorithm under the same conditions are used as the baseline case for comparison. For additional comparisons, a fuel consumption measure is computed as the summation of the absolute commanded actuations (i.e., ftzyxiipdtuF0,,)( for the position and ftiiadtTF02,1)(9180.0JKNFp/0815.0 for the attitude). For the baseline case, the cost functional value is and the fuel consumed by the position and attitude control maneuvers are and respectively. Kg MNFa7185.9 The genetic algorithm parameters used in the simulation are listed in Table 4-4. The parameters range in SMC + BLC is defined by experience (e.g., the initial relative position and velocity error, and attitude error et al). k is related to thruster bandwidth, which represents how faster the thruster can force the performance switch between the sliding surfaces, is the corresponding proportional gain, and determines how smooth the performance will be. [10;10;200][500;500;1000]pk [0;0;0][1;1;1]p [10;1;0][100;60;20]ak

PAGE 150

132 [1;1;0][8;7;3]a [0.001;0.001;0.001][0.1;0.1;0.1]p [0.01;0.001;0.001][0.1;0.01;0.1]a Table 4-4 GA parameters definition and the data used Parameter Physical Meaning Parameter Physical Meaning 18_Paran Number of Parameters need to be optimized 32_ Genn Number of Generation 40bl Bit Length 9.0_ Cxp Probability of Crossover 5_Pops Population Size 05.0_ Mutp Probability of mutation Under the same simulation condition (GA parameters, initial orbits, control parameters, and thruster maximum levels), the simulation is repeated for nine times. The results are listed in Table 4-5. The value of the cost functional for all nine cases is smaller than the baseline case in which the controller parameters were somewhat arbitrarily chosen. The feature that all nine cases give different cost value is due to the fact that local minimum are found for each case. The minimum cost functional in this table (case 2) 0.2505 is reduced by 71.8% compared with the baseline case (in the sense of difference percentage). Furthermore, all the cost values in Table 4-5 are on the same order, which implies that many minimum or sub minimum points exist in the chosen range. The fuel consumption in relative position and attitude are (reduced by 62.4%) and KgKNFp/0502.0 MNFa 48256.3 (reduced by 64.1%). Case 2, which showed the most improvement in the cost functional, was chosen from Table 2 and used to simulate the systems response and the results for the relative position are shown in Figures 4-58 through 4-63. As shown in Chapter 2, the z-direction dominates the relative position error, and since the position errors were equally weighted,

PAGE 151

133 the best improvements were observed in the z-direction (Figure 4-62). The control forces in all three directions show the same trends, where the genetic algorithm results are better in the z-direction and are not necessarily better in the xand y-directions. Table 4-5 Simulation results of GA based SMC + BLC Cases pk p ak a p a Cost 1 [149.8754; 156.4414; 920.8560] [0.3913; 0.4234; 0.1958] [19.7441; 25.4264; 7.4341] [3.7798; 1.8269; 2.1746] [0.05; 0.0497; 0.0190] [0.0530; 0.0064; 0.0629] J = 0.3870 2 [50.9523; 64.3216; 979.7636] [0.2378; 0.902; 0.3280] [33.8714; 44.6930; 2.6078] [2.3105; 1.0567; 1.4913] [0.0583; 0.0057; 0.0046] [0.0752; 0.0032; 0.0990] J = 0.2505 3 [459.8061; 348.2180; 989.0553] [0.2629; 0.9223; 0.1755] [70.8052; 32.7032; 14.6047] [1.1438; 1.0650; 1.9601] [0.0190; 0.0693; 0.0061] [0.0633; 0.0084; 0.0607] J = 0.2755 4 [442.3733; 99.9873; 797.0193] [0.5562; 0.5722; 0.2591] [52.5615; 40.0959; 5.1101] [2.4270; 2.5129; 1.7649] [0.0024; 0.0407; 0.0067] [0.0264; 0.0064; 0.0434] J = 0.3264 5 [126.2625; 381.5813; 954.1601] [0.5351; 0.2559; 0.6135] [96.6122; 28.8326; 4.3939] [4.7652; 1.0547; 2.7337] [0.00867; 0.0770; 0.0075] [0.0514; 0.0064; 0.0719] J = 0.3340 6 [472.0258; 135.2275; 965.5037] [0.3823; 0.4286; 0.0918] [86.6887; 43.5448; 4.5561] [6.3803; 1.0303; 1.6382] [0.0454; 0.0623; 0.0146] [0.0879; 0.0075; 0.0359] J = 0.3002 7 [208.9834; 337.9914; 926.3969] [0.6837; 0.4637; 0.5561] [36.7344; 49.6019; 8.8869] [5.6633; 1.4342; 2.0515] [0.0162; 0.0364; 0.0145] [0.0212; 0.0095; 0.0860] J = 0.3149 8 [104.5852; 94.1560; 934.3774] [0.2007; 0.7149; 0.1404] [66.6932; 58.5784; 9.4989] [2.7227; 1.8357; 2.0586] [0.0823; 0.0017; 0.0037] [0.0143; 0.0029; 0.0461] J =0.2564 9 [498.2046; 32.2731; 772.4431] [0.0024; 0.1199; 0.3315] [25.1473; 56.2396; 13.6366] [6.4325; 1.6175; 0.9156] [0.0609; 0.0063; 0.0898] [0.0836; 0.0760; 0.0050] J =0.2661

PAGE 152

134 Figure 4-58 The relative position error in xdirection. Figure 4-59 Control force in x-direction. Figure 4-60 The relative position error in ydirection. Figure 4-61 Control force in y-direction Figure 4-62 The relative position error in zdirection. Figure 4-63 Control force in z-direction Figures 4-64 through 4-67 show the simulation results for the attitude errors. The same reason (Azimuth error angle dominate the relative attitude error) as the relative position simulations shown above, for the settling time of the case with GA in the elevation error angle is much longer than that of the case without GA, the required

PAGE 153

135 torques are less. For the control torque commands shown in Figure 4-67, the different saturation level is coming from the different boundary layer parameter a Figure 4-64 Elevation error angle. Figure 4-65 Torque in body axis 1. Figure 4-66 Azimuth error angle. Figure 4-67 Torque in body axis 2. Up to now, the low level controller has been designed. The sliding mode controller (SMC) was chosen as the control methodology used in this dissertation for the attitude and orbit dynamics based on the discussions of chapter 3. However, the SMC demonstrated some inappropriate responses (i.e., chatter and excessive fuel consumption) which were addressed in this chapter. A boundary layer control (BLC) was added to BLC to mitigate the chattering phenomena, and then an optimization methodology was used to select the best controller parameter. Because the system is highly nonlinear, a genetic algorithm was applied for the optimization and the fuel consumption is improved more than 60%.

PAGE 154

CHAPTER 5 RECONFIGURATION CONTROLLER IN FORMATION FLYING 5.1 Introduction There are many advantages of the formation flying system, such as cost reduction, increasing baseline observation, reliability, and mission flexibility. In order to achieve these goals, especially reliability and mission flexibility, multiple vehicles maintenance and reconfiguration control is a key technology that is required. Within the multiple vehicles control issues, the formation architecture is the first topic to be considered. In the literatures, there are four types of architectures; they are shown in Figure 5.1 through Figure 5.4 [84]. Figure 5.1 shows the absolute architecture in which each satellite communicates with the ground station. The ground station is in charge of the decision-making. Due to the distance between the Earth and the satellite around it, using this feedback strategy to deal with noise and errors is not acceptable, because the feedback delays associated with communication. In order to avoid the feedback, every satellite needs to have a long redundancy data frame to check and correct the error bits. This will reduce the data throughput and may actually degrade the system stability. Seen in Figure 5.2, the leader-follower architecture, the Leader satellite communicates with the ground base, whereas all the other satellites (follower satellites) only exchange the data information with the leader. Because the communication between the leader and the ground base occurs less often, the data rate of the system is increased and then the control precision may be improved. 136

PAGE 155

137 Figure 5.1 Absolute Figure 5.2 Leader-follower Structure Figure 5.3 Absolute and interconnection Figure 5.4 Sub formation The absolute and interconnection architecture (shown in Figure 5.3) is a more complicated version of the absolute one. Every satellite in the formation exchanges information with each other and also with the ground station. In this architecture, any one can be the leader depends on the mission and fault conditions. The sub formation architecture shown in Figure 5.4 deals with the system that has many satellites, which makes the system to be easier to manage under a hierarchical structure. The top satellite is the leader, which communicates with the ground-station and sub-leaders. The sub-leader exchanges information with the followers and the leader. The architectures in Figure 5.3 and Figure 5.4 are complicated in the software and hardwares design and implementation although they may have more robustness and better flexibility.

PAGE 156

138 Considering the trade off among the flexibility, complexity, and data throughput, the leader follower architecture is chosen for this dissertation. As an example, the architecture of a three-satellite in-plane formation is discussed in this section. As seen in Figure 5.5, the basic formation is composed of 3 satellites flying in the same orbit plane, one is the leader (L) and the other two are followers (F1 and F2). This plane is called the mission plane here. All the satellites in the mission plane are controlled to point to the same point, i.e., the nadir point of the leader The leader satellite is assumed in perfect working condition during the whole mission. In the case, when the follower satellite fails, two redundant satellites (R1 and R2), which are initially flying in an orbit slightly inclined to the mission plane, will take place of F1 and (or) F2. This inclined plane is referred to as the redundant plane here. LN As an example, in Table 5-1, the initial orbit information for each satellite of this formation configuration is listed. The leader and two followers are flying in the same equatorial plane, which is the main plane, with 0.1 degrees separation between them. Two redundant satellites are flying in a slightly inclined orbit with 0.1 degrees inclination, which has true longitude. Because all the satellites have the same orbit radius and different true longitude at epoch, two redundant satellites are not crossing the mission plane at the same time as F1 or F2 in the open loop. Furthermore, in Table 5-1, the follower and redundant satellites initial relative distances from the leader are also listed. The small inclination and true longitude are chosen in order to get an order of 10 km in the relative position. 00.05

PAGE 157

139 Figure 5.5 3-satellite formation architecture with two redundant satellites Table 5. 1 The initial orbit information of the formation satellites Leader Follower 1 Follower 2 Redundant 1 Redundant 2 Right Ascension 0o 0o 0o 0o 0o Inclination 0oi 0oi 0oi 0.1oi 0.1oi Eccentricity 0e 0e 0e 0e 0e Semi-major 6878ak mmmmm 6878ak 6878ak 6878ak 6878ak Initial true anomaly 270of 269.9of 270.1of 269.95of 270.05of Argument of periapsis 90o 90o 90o 90o 90o True longitude 0os 0.1os 0.1os 0.05os 0.05os Initial Distance to the Leader 12.00437307157397 km 12.00437307157397km 6.00217796530315 km 6.00217796530315 km In this Chapter, first, some assumptions and considerations are listed for the safe hold and fuel consumption. After that, the perceptive frame theory is introduced in Section 5.3. In Section 5.4, the corresponding event-driven mapping is derived and the low level controllers based on SMC and BLC are designed. Finally, the simulation results based on perceptive frame, SMC and BLC are shown. A three-satellite formation

PAGE 158

140 mentioned above is used as an example to demonstrate the formation maintenance and reconfiguration. 5.2 Formation Reconfiguration Strategy Because of the mission changes and error conditions, the formation of the system may need to be changed. When Follower 1 or Follower 2 fails in the mission, the threesatellite configuration is changed. The fo llowing assumptions are made here for the formation reconfiguration strategy. Assumption 1: The leader satellite works perfect and never fails in the simulation. Assumption 2: When a follower satellite works fine except that it cannot communicate with the leader satellite properly, this follower satellite automatically drives itself to a safe orbit, where it does no t collide with all the other satellites. Assumption 3: if the follower satellites thruster does not work, there exist some other methods (for example, a satellite tow truck), which can be used to move this defective satellite to a safe orbit. Assumption 4: In order to minimize the fuel consumption, all the switches happen when the redundant satellite is passing thr ough the ascending or descending node of the redundant plane depends on the location of the failure occurs. Assumption 5: In order to minimize the fuel consumption, redundant satellites are in a free flying status before the switching for reconfiguration. This means, when the initial three satellites formation works fine the redundant satellites are in an unforced condition. Only when a reconfiguration is commanded will the redundant satellites be controlled. Furthermore, because the satellite s in the mission planes satellites are all equatorial and circular with the same ra dius, fuel consumptions only exist in the reconfiguration transient stages and orbit altitude maintenance.

PAGE 159

141 In the case when the reconfiguration of the formation is needed, the nearest redundant satellite is driven by the controller to the new position. For example, in Figure 5.6 (and Figure 5.7), when the follower 2 (1) fails, the redundant satellite 2 (1) is driven to the main orbit plane and has a true longitude angle less than the previous follower satellite 2 (1) for the safe consideration. The orbital information of the reconfigured redundant satellites, including the desired states, is listed in Table 5-2. Figure 5.6 If F 2 fails

PAGE 160

142 Figure 5.7 If F 1 fails Table 5-2 Orbit information of R1 and R2 in figures 3 and 4 Redundant 1 Redundant 2 Right Ascension 0o 0o Inclination 0oi 0oi Eccentricity 0e 0e Semi-major 6878ak m 6878ak m True anomaly 269.95of 270.05of Argument of periapsis 90o 90o True longitude 0.05os 0.05os Initial Distance to the Leader 6.00218710715283 km 6.00218710715284 km Desired states (Relative position, relative velocity, angle, angular velocity) [-0.00261894810046; -6.00218653578698; 0;0;0;0;0;0;0;0; -1;0] [-0.00261894810046; 6.00218653578699; 0;0; 0;0;0;0;0;0; -1;0] Desired absolute position [6877.99738105190; -6.00218653579;0]km [6877.99738105190; 6.00218653579;0] km When both F1 and F2 fail, R1 and R2 replace F1 and F2. The orbit information for R1 and R2 are the same as Table 5-2. The reconfiguration switching time is at ascending or descending node for R1 and R2 separately.

PAGE 161

143 5.3 Perceptive Frame Control (PFC) Based on the previous sections, formation flying problem can be separated as centralized or decentralized control problem. The centralized controller uses the absolute communication architecture. Although the control performance may be better than the decentralized one, the communication load is heavy. The reconfiguration is not flexible and smooth. On the other hand, for the leader-follower or decentralized architecture, the leader determines when and where the reconfiguring maneuver should be. In general, the leader-follower controller only deals with the low level continuous controller (e.g., maintenance and tracking), and does not reflect the high level events. This is because, in most cases, the controllers independent variable is time, which is not directly or easily related to the event. Perceptive frame (PF) is the connection between the high level discrete event and the low level continuous controller [38, 39]. Therefore, individual vehicles in the formation can have different control strategies. The whole formation is coordinated by the high level discrete event. Then, the time-dependant feedback control becomes event-driven, or formation dependant controller. So, it is easy to implement in the distributed system, and changing the number of the satellites in the formation will not affect most parts of the system. Two basic concepts in the PF are the action reference and reference projection. The action reference is the key parameter, which represents the configuration and reconfiguration of the mission in a most convenience way. Different mission or task of the formation has different action reference. As an example in references [39, 85], in a 2-dimension, in plane and circular formation flight system, where satellites are separated by an angle the most convenience action reference is the phase angle This angle

PAGE 162

144 represents satellites position in the circular orbit. Therefore, the initial action reference for the leader is 1L and that of the follower is 1L 1 where 1 is the angle separation between the leader and the follower. When the formation is changed, the leaders action reference becomes 2L and that of the follower is 2 2 L Not like the time dependant parameter, the action reference is determined by the mission level commands. The reconfiguration trigger is an event instead of time. )t 1ij dLy (,)x ()) s s The reference projection is a mapping from the action reference to the controllers independent variable, such as time. The projection can be continuous and discontinuous, depends on particular problems. Just as the above example, a mapping from the action reference to the controllers time is t where is the angular velocity. Therefore, controller can be written as uu (uu ( ). Under this reference projection, the time dependant controller becomes the event dependant one. A perceptive frame has four steps in the formation controllers design. The first step is to generate the action reference and the desired path for every vehicle i in the formation. The desired path needs to be a function of the action reference as where is the leaders desired trajectory, and s( diy ()()()PdidLjjysyses()ijs )s ()jes are the unit vector and each vehicles coordinate in the leaders frame separately. Step 2 designs the controller uu t for each vehicle in the formation separately. Step 3 defines the reference projection (s ) x which satisfy ( d x This reference projection maps the desired states to the action reference. Furthermore, constructing mapping function from the action reference to the independent variable t in the controller

PAGE 163

145 ()svt For the simple case, the function vt can be continuous. But in most cases, the reference projection function is discontinuous or discrete. Step 4 involves the non-time based control low from the reference projection as Figure 5.8 shows the basic structure for the perceptive frame controller. The stability theorem is discussed and proven in reference [72]. () 1(,)(,(()))uxtuxvx Figure 5.8 Perceptive frame controllers structure 5.4 Controller Design In this section, an example is given to validate that the perceptive frame and genetic algorithm based sliding mode control is capable of maintaining and re-configure the formation flight system. First the perceptive application is addressed for the particular formation mission. Then, the formation maintenance control (FMC) and formation reconfiguration control (FRC) are designed. 5.4.1 Perceptive Frame Controller For the mission mentioned in Sections 5.2 and 5.3, both the relative position and attitude are controlled. The mission decision provides the desired configuration code based on the ground commands and error handling information. SMC & BLC controller

PAGE 164

146 maintenances and reconfigures the formation in the low level continuous sense. The perceptive frame constructs the connection between the high-level configuration code and the low-level SMC & BLC controller. Because the relative attitude controller is coupled with the real time relative position, then any time, the relative positions maneuver triggers the attitude controller to give a corresponding torque commands. Therefore, the most convenience action reference is the relative position between the leader and follower. However, based on Chapters 2, 3 and 4, it is very difficult to find the analytical relation between the relative position and controllers independent variable, such as time. Therefore, in this Chapter, an indirect event, true anomaly is chosen as the action reference. Also, the reference projection is a discrete function as in Eq. (5.1). equals zero represents the case that the control is rest. equals one and two are the cases the controller been triggered at ascending and descending nodes respectively. Lf Flag Flag f is the true anomaly of the satellite. 1,2,3R represents the different fault cases and 0R is the normal case. 012Flag when if 01,2,31,2,3RRR and (5.1) 0180180360oooff 5.4.2 Formation Maintenance Controller The Formation Maintenance Controller (FMC) is designed, which tries to maintain the desired states. The parameters for the sliding mode and boundary layer controllers are selected based on the genetic algorithm simulation. The chosen parameters are shown in Table 5-3. Here, no control effort is acting on the redundant satellites, because, in order to minimize the fuel consumption, in the maintenance phase, all the redundant satellites are in free flying status.

PAGE 165

147 Table 5-3 Formation maintenance controllers for follower 1 and follower 2 FMC parameters Follower 1 Follower 2 p K [50.6176; 234.7942; 818.5500]T [326.9580; 320.5726; 417.1905]T p [0.7698; 0.0782; 0.3595]T [0.4065; 0.7695; 0.6331]T a K [16.2940; 52.0524; 5.2412]T [32.2423; 50.9221; 14.3360]T a [4.8082; 1.0123; 0.0164]T [5.8961; 1.0139; 1.0360]T p [0.0053; 0.0530; 0.0665]T [0.0194; 0.0170; 0.0958]T a [0.0712; 0.0011; 0.0916]T [0.0568; 0.0011; 0.0596]T 5.4.3 Formation Reconfiguration Controller When F1, F2, or both have failed, the formation is reconfigured as shown in Figure 5.6 and Figure 5.7. Different control parameters based on the genetic algorithm are selected and listed in Table 5-4. Because the redundant satellites are driven to replace the follower satellites, these controllers are named as Formation Reconfiguration Controller (FRC). Table 5-4 Formation reconfiguration controllers for redundant 1 and redundant 2 FRC parameters Redundant 1 Redundant 2 p K [446.5047; 12.5296; 405.3798]T [432.6077; 132.2517; 921.2565]T p [0.0662; 0.7089; 0.3956]T [0.6964; 0.2042; 0.3406]T a K [36.1409; 57.0358; 15.3765]T [20.3617; 49.0118; 0.3756]T a [2.1098; 1.4490; 2.7421]T [7.9901; 1.3662; 2.1134 ]T p [0.0854; 0.0046; 0.0077]T [0.0769; 0.0539; 0.0593]T a [0.0906; 0.0046; 0.0709]T [0.0198; 0.0024; 0.0957]T 5.4.4 SIMULINK Block Diagram Figure 5.9 gives the SIMULINK diagram for the formation reconfiguration controller. The clock block is used to broadcast the leaders true anomaly. The configuration code is passed from the decision block in the satellites. The decision block generates the configuration code based on information coming from the communication

PAGE 166

148 and sensor systems. Nevertheless, here, the configuration code is simulated using MATLAB. Table 5-5 listed the simulated configuration code used in the simulation. After the configuration code passed to the Global Variable block, corresponding global control parameters and desired states are broadcasted. Only related models within four satellites block update their control parameters and desired states. For example, if the configuration code received is 1, which represents the case when the follower 1 failed, the redundant satellite 1s control parameters and desired states are broadcasted. Because only the leader / redundant 1 block has the global variables related to these control parameters and desired states, only this block updates its information. Based on the controllers parameters and desired states, the leader / follower or leader / redundant block calculated the closed loop response and send them to the respective output block. The output includes error angles, relative position, force and torque commands, and leaders absolute position. Figure 5.9 SIMULINK program diagram

PAGE 167

149 Table 5-5 Configuration code for in-plane formation flight simulation Configuration Code Formation Flight Diagram 0 All the follower satellites work (Figure 5.5) 1 Follower 1 failed 2 Follower 2 failed 3 Both followers failed 5.5 Results Total three cases are simulated, which ar e coded as 0-1-3, 0-2-3 and 0-3. For example, 0-3 represents follower 1 and 2 of the formation system both fail. Here, only 03 simulation results are shown. The relative positions are shown in Figure 5.10 through Figure 5.12. In all these relative position difference results, the relative position for the followers 1 and 2 are constants due to the fact that the formation is in plane, the followers one and two are initially in the desired states and the orbit is circular. Redund ant satellites 1 and 2 are free flying in the beginning as shown in Figur e 5.10 (a) through Figure 5.12(a). When the satellites reach the ascending node separately, the formation is reconfigured due to the simulated error conditions. At this time, in local directions, R1 and R2 track the new desired states corresponding to the configur ation commands. As shown in Figure 5.10 (b), Figure 5.11 (b) and Figure 5.12 (b), the switching happens when the redundant satellites cross the mission orbit, and R2 crosse s the mission orbit first. Just as shown in Figure 5. 11 (a), the stable states of redundant satellites are not th e same as the follower due to the different true longitude condition.

PAGE 168

150 Figure 5.10a Figure 5.10b Figure 5.10 Relative position in x local frame. (a) whole range, and (b) from 5670 7500 second. Figure 5.11a Figure 5.11b Figure 5.11 Relative position in y local frame, (a) whole range, and (b) from 5500-10000 seconds.

PAGE 169

151 Figure 5.12a Figure 5.12b Figure 5.12 Relative position differences in z local frame, (a) 0-800 seconds, and (b) 5650-6500 seconds The average settling time is calculated by 11nSettlingiSettlingitn t, where n is the number of the transient stages and t represents the settling time for state response. For example, seen in Figure 5.10, the settling times of the relative position in iSettling thi x for the redundant satellites are approximately 1322 and 723 seconds (Figure 5.10 (b)). In Figure 5.11 (b), the settling times in the y direction are both 1500 seconds. In Figure 5.12(b), the settling times are 173 and 422 seconds. Therefore, for the relative position, the settling time is approximately 940 seconds, which is nearly 17% of the period.

PAGE 170

152 Figure 5.13a Figure 5.13b Figure 5.13 Force commands in x local frame. (a) whole range, and (b) 5600-10000 seconds Figure 5.14a Figure 5.14b Figure 5.14 Force commands in y local frame. (a) whole range, and (b) 5600-9000 seconds

PAGE 171

153 Figure 5.15a Figure 5.15b Figure 5.15 Force commands in z local frame. (a) 0-800 seconds, and (b) 5900-6300 seconds Figure 5.13 to Figure 5.15 show the force commands. Because the formation is in a circular orbit, there is no control in the follower satellites. Using the sliding mode and boundary layer controllers, the control commands for the redundant satellites are small. Just as the relative position response, the average settling time is large. In Figure 5.15(b), the boundary layer effects of redundant 1 can be seen. First, the state is out of the boundary layer; therefore, the saturation effect is shown. When the state is within the boundary, smooth response is obtained. Here, two redundant satellites have different saturation level due to the different boundary layer parameters are used here. In Table 5-4, pz for redundant 1 is 0.0077, whereas that of redundant 2 is 0.0593. This is the reason, why redundant 1 goes to a lower saturation level than redundant 2. The total fuel consumption for the relative position control is calculated by 211mnijpositionijfuelF where [1,2,...,]i m represents the satellites, and here, 1,2i is the follower 1 and 2 and 3,4i represent the redundant satellites 1 and 2.

PAGE 172

154 [1,2,...,]j n denotes the commands index. [,,]TijijxijyijzFFFF is the th j force command for i satellite. For the relative position control, the fuel consumption is th -5=3.33157287423550210 KN/Kg Positionfuel Figure 5.16a Figure 5.16b Figure 5.16 Elevation error angle. (a) whole range, and (b) 5670-15000 seconds Figure 5.17a Figure 5.17b Figure 5.17 Azimuth error angle. (a) whole range, and (b) 5900-7000 seconds From Figure 5.16 through Figure 5.17, the error angles time responses are shown. The azimuth error angle is bounded, and has a settling time of 523 seconds, whereas, the time response for the elevation error angle is very slow and has a high frequency

PAGE 173

155 oscillation. The average settling time for the elevation error angles is large than one period (5677 seconds). The control resolution for the elevator error angles is on the order of radians, whereas, that of the azimuth error angle is on the order of O. The elevation error angle has stable states around zero. But the stable states of the azimuth error angles for R1 and R2 are around 6(10)O 4(10) 47.6410 3 and radians, which is less than that of F1 and F2 ( 48.110 1.0910 radians). This is because the true longitude of the redundant satellites is less than that of the follower. The torque commands for the attitude model are shown in Figure 5.18 Figure 5.20. The fuel consumption is 0.04228326917338NMAttitudefuel The torques for F1 and F2 in body axis 1 are zeros, because the followers are already in the main orbit, Whereas the torques for F1 and F2 in body axis 2 are driven to zeros from initial states. Also, just like the relative position control effort, the saturation levels for R1 and R2 are different. Figure 5.18 Torque commands in body axis 1. Figure 5.19 Torque commands in body axis 2.

PAGE 174

156 Figure 5.20 Torque commands in body axis 2, from 5670 7000 seconds 5.6 Summary In this Chapter, the perceptive frame concept is used to construct the relation between the high level discrete event and the low level continuous controller. A three-satellite formations configuration and reconfiguration is discussed.

PAGE 175

CHAPTER 6 CONCLUSION AND FUTURE WORK 6.1 Conclusion In this research, the dynamics and control for the satellite formation flying system are proposed based on the leader-follower structure. First, a brief introduction of the satellite formation flying system is discussed, including formation flying systems properties, potential advantages, challenges, and the focused research scopes in this dissertation. Second, due to the lack of closed loop simulation environment in STK and Free Flyer, the relative motion dynamics model is developed. In order to have a model which can be used in any closed orbit cases, including circular and elliptic reference orbits, the true anomaly of the leader satellite is chosen as the independent variable instead of time. Based on the leader follower structure, the multi-agents dynamics system is simplified as a two-body problem. Different to most of the literatures to date, coupled relative attitude dynamics model is derived related to the real time relative position. Due to the mission requirement that all the satellites in the formation should have a pointing to the nadir of the leaders orbit managed collectively, the desired attitude is coupled with the real time position. Based on the spherical triangle theory, the error angles between the actual and desired states are calculated. Simplifying the derived generalized dynamics model, the often appeared models, such as C-W model, is introduced as an special cases. Atmospheric drag and J 2 perturbation from the Earth Oblateness are also introduced into the dynamics model. The effects of the drag and perturbation are studied. 157

PAGE 176

158 The formation flying control system is a hybrid system, with a low-level continuous system and a high-level discrete system or a low-level configuration maintaining system and a high-level reconfiguration system. Therefore, in Chapter 3, the low level continuous systems control problem is studied. Linear controllers, such as LQR and /H and nonlinear controllers, such as feedback linearization and sliding mode control, are designed based on the C-W equation and fully nonlinear model, respectively. Simulation results are compared, and finally, sliding mode control strategy is chosen as the base method for both orbit and orientation control in this dissertation. Furthermore, in order to reduce the chattering phenomena, fuel consumption and improve the performance, boundary layer control and genetic algorithm are used and combined with the sliding mode control. As a supplement issue, the discrete PID controller is designed for the purpose of orbit altitude degradation mitigation due to the atmospheric drag and J 2 perturbation. In the end, the high-level reconfiguration control is discussed in Chapter 4. In this Chapter, the orbit mission is designed with the consideration of fuel consumption, collision avoidance, fault handling, and simplicity for implementation. Configuration maintaining and reconfiguration controllers are designed based on sliding mode control, boundary layer control and genetic algorithm. 6.2 Future Work The research topics for the satellite formation flying system is open ended. Some areas, which will benefit the formation flying system, are listed below. 1. Precision control of the relative position and attitude to the level of less than micro-meter and micro-arc-second. In this dissertation, combined sliding mode control and boundary layer control is capable of controlling the system to this level of precision. However, all the actuators / thrusters have the assumption that they can support infinite

PAGE 177

159 resolution. But in reality, this kind of thruster cannot be provided. Therefore, how to design the precisely controller under the limited thruster resolution is one of the future research focus. 2. Orbit design. In space mission, the fuel cons umption is a very important issue. How to design the orbits for the formation sy stem for the purpose of reducing the fuel consumption, avoiding collisions, and making the launching as simple as possible are also a key technology, which need to be solved in the near future. 3. Hybrid system. In this dissertation, the basic idea of the perceptive frame is introduced in the configuration and reconfi guration control. However, the theoretical analysis of the perceptive frame has not been proven to date. 4. Safe hold algorithm. The simple safe hold id ea has been discussed in this dissertation for a 3-satellite formation. In the future, a fu ll-concept safe hold algor ithm is need to be proposed, including health detection, error handling algorithm and corresponding control actions. 5. Mini-satellites hardware manufact ure. This includi ng the light-weight communication, thruster, and sensor hardware.

PAGE 178

APPENDIX QUATERNION AND LIE DERIVATIVE A.1 Quaternion From [Ref. 56], rotation about an arbitrary axis fixed in body frame is defined to be quaternion parameters: 11sin/2qe 22sin/2qe 33sin/2qe 4cos/2q and 123[,,]Tqqqq Where the arbitrary axis is: [,. Rotation angle is: 123,]Teee These four parameters are not independent of each other. The magnitude of this vector is 1: 241Tqqq The relation between direction cosine matrix and quaternion is: 22231234132422123413231422132423141224412()2()2()2()12()2()2()2()12()()22tTqqqqqqqqqqCqqqqqqqqqqqqqqqqqqqqqqqIqqqQ Where 323121000qqQqqqq 160

PAGE 179

161 And 4112233(1)/4qCCC 233231134122114CCqCCqCC if q 40 A.2 Lie Derivative [Ref 56] Given a smooth scalar function the gradient of h is: ()hx hh x Also let to be the vector function mapping f n R to n R then Jacobian of f is: f f x and () iijj f f x Definition B.2.1: Lie Derivative Let hR be a smooth scalar function, and :nR :n n f R R be a smooth vector field on n R then the Lie derivative of with respect to h f is a scalar function defined by Thus, the Lie derivative is simply the directional derivative of along the direction of the vector fLh hf fLh h f

PAGE 180

LIST OF REFERENCES [1]. Steve Garber, NASA History Web Curator, February 21, 2003 (last update), http://www.hq.nasa.gov/office/pao/History/sputnik October 28, 2003 (last access). [2]. Kate R.Hartman, Cheryl J. Gramling, Taesul Lee, David A. Kelbel, and Anne C. Long, Relative navigation for spacecraft formation flying, AAS/GSFC 13 th International Symposium on Space Flight Mechanics, AAS 98-379, Greenbelt MD, May 11-15, 1998. [3]. P.K.C.Wang, J.Yee, and F.Y.Hadaegh, Synchronized rotation of multiple autonomous spacecraft with rule-based controls: experimental study, Journal of Guidance, Control, and Dynamics, Vol. 24 No. 2, pp. 352-360, March-April 2001 [4]. Christoper Kitts, and Freddy Pranajaya, Emerald: an experimental mission in robust distributed space systems, 13 th AIAA/USU Conference on Small Satellites, SSC-99-VI-5, Logan, Utah, August 23-26, 1999. [5]. NASA New Millennium Program, September 17, 2003, (last update), http://nmp.jpl.nasa.gov/ October 28, 2003 (last access). [6]. Andrew W. Lewin, Creating large space platforms from small satellites, 13 th AIAA/USU Conference on Small Satellites, SSC-99-VI-6 Logan, Utah, August 23-26, 1999. [7]. Hans F. Meissinger, John Collins, Gwynne Gurevich, and Simon Dawson, Low cost, minimum size satellites for demonstration of formation flying modes at small, kilometer-size distances, 13 th AIAA/USU Conference on Small Satellites, SSC-99-VI-3, Logan, Utah, August 23-26, 1999. [8]. Jeff Ward, Susan Jason, and Martin Sweeting, Microsatellite constellation for disaster monitoring, 13 th AIAA/USU Conference on Small Satellites, SSC-99-V-2, Logan, Utah, August 23-26, 1999. [9]. Richard Blomquist, Solar blade nano-satellite development: heliogyro deployment, dynamics and control, 13 th Annual AIAA/USU Conference on Small Satellite, SSC-99-VII-2, Logan, Utah, August 23-26, 1999. [10]. C.D.Rayburn, H.E.Spence, H.E.Petschek, M.Bellino, J.Vickers, and M.Murphy, Constellation pathfinder: a university nanosatellite, 13 th AIAA/USU Conference on Small Satellites, SSC-99-V-1, Logan, Utah, August 23-26, 1999. 162

PAGE 181

163 [11]. Dave Speer, Phyllis Hestnes, Mark Perry, Bill Stabnow, The new millennium program EO-1 mission and spacecraft design concept, IEEE Aerospace Conference pp. 207-227, Vol. 4, Snowmass, CO., February 1-8, 1997. [12]. S. Horan, B. Anderson, B. Underhill, A. Friedman, J. Wong, H. Reed, E. Hansen, and D. Rodier, Three corner sat constellation --New Mexico State University: communications, LEO telecommunications services, intersatellite communications, and ground stations and network, 13 th Annual AIAA/USU Conference on Small Satellite, SSC-99-VI-7, Logan, Utah, August 23-26, 1999. [13]. Frederic Teston, Richard Creasey, J. bermyn, D. Bernaerts, K. Mellab, PROBA: ESAs autonomy and technology demonstration mission, 13 th Annual AIAA/USU Conference on Small Satellite, SSC-99-V-8, Logan, Utah, August 23-26, 1999. [14]. Peter Rossoni, and Peter V. Panetta, Developments in nano-satellite structural subsystem design at NASA-GSFC, 13 th AIAA/USU Conference on Small Satellites, SSC-99-V-3, Logan, Utah, August 23-26, 1999. [15]. Harold Dahl, Wall Eliuk, Glen Rumbold, and Randolph Shelly, ACE-a Canadian small science satellite mission, 13th AIAA/USU Conference on Small Satellites, SSC99-V-7, Logan, Utah, August 23-26, 1999. [16]. David A. Weidow, John Bristow, NASA/DoD university nano-satellites for distributed spacecraft control, 13th Annual AIAA/USU Conference on Small Satellite, Logan Utah, August23-26, 1999. [17]. Ruth Moser, David Collins, Alok Das, Robert Ferber, George Jaivin, Richard Madison, Joseph Smith, and Michael Stallard, Novel missions for next gerneration microsatellites: the results of a joint AFRL-JPL study, 13 th AIAA/USU Conference on Small Satellites, SSC-99-VII-1, Logan, Utah, August 23-26, 1999. [18]. Zsolt Kiraly, Brian Engberg, Franz Busse, and Robert J. Twiggs, The ORION microsatellite: a demonstration of formation flying in orbit, 13 th Annual AIAA/USU Conference on Small Satellite, SSC-99-VI-8, Logan, Utah, August 23-26, 1999. [19]. Blaise Morton, and Nicholas Weininger, Collective management of satellite clusters, AIAA Guidance, Navigation and Control Conference, pp. 1576-1584, Portland, Oregon, August 1999. [20]. Manop Aorpimai, Phil Palmer, Alex da Silva Curel, and Surrey space Centre, Phase acquisition and formation keeping of a new power consumption monitoring satellite constellation, 13 th Annual AIAA/USU Conference on Small Satellite, SSC-99-VI-2, Logan, Utah, August 23-26, 1999.

PAGE 182

164 [21]. Richard H. Vassar, and Richard B.Sherwood, Formation keeping for a pair of satellites in a circular orbit, Journal of Guidance, Control, and Dynamics, Vol. 8, No. 2, pp. 235-242, March-April 1985. [22]. David C. Redding, and Neil J. Adams, Linear-quadratic station keeping for the STS orbiter, Journal of Guidance, Control, and Dynamics, Vol. 12, No. 2, pp. 248-255, March-April 1989. [23]. R.J.Sedwick, D.W.Miller, and E.M.C.Kong, Mitigation of differential perturbations in formation flying satellite clusters, The Journal of the Astronautical Sciences, Vol. 47, No. 3 and 4, pp.309-331, July-December 1999. [24]. William Singhose, Tarunraj Singh, and Warren Seering, On-off control of flexible spacecraft with specified fuel usage, American Control Conference,pp. 2308-2312, Vol. 4, Albuquerque, New Mexico, June 4-6, 1997 [25]. Hui Yan, and Hongxin Wu, Optimal low-thrust earth-moon targeting strategy for n-body problem, Journal of Guidance, Control, and Dynamics, Vol. 24 No. 3, pp. 626-628, May-June 2001 [26]. Qiguo Yan, Vikram Kapila, and Andrew G. Sparks, Pulse-based periodic control for spacecraft formation flying, American Control Conference, pp. 374-378, Vol. 1, Chicago, IL, June 28-30, 2000. [27]. R. K. Yedavalli; A. G. Sparks, Satellite formation flying control design based on hybrid control system stability analysis, American Control Conference, pp. 2210-2214, Vol. 3, Chicago, IL, June 28-30, 2000. [28]. S. M. Veres, S. B ,Gabriel, E. Rogers, D. Q. Mayne, Analysis of formation flying control of a pair of nano-satellites, Decision and Control, pp. 1095-1100, Vol. 2, Orlando, FL, December 2001. [29]. S. Schweighart, R. Sedwick, A perturbative analysis of geopotential disturbances for satellite cluster formation flying, IEEE Aerospace Conference, pp. 1001-1019, Vol. 2, Arlington, Virginia, June 25-27, 2001. [30]. S. R. Starin, R. K. Yedavalli, Sparks, A.G., Design of a LQR controller of reduced inputs for multiple spacecraft formation flying, IEEE American Control Conference, pp. 1327-1332, Vol 2, Arlington, Virginia, June 25-27, 2001. [31]. De Queiroz, V.Kapila, and Q.Yan, Adaptive nonlinear control of satellite formation flying, AIAA Guidance, Navigation and Control Conference, pp. 1596-1604, Portland, Oregon, August 1999. [32]. de Queirox, M.S.; Yan, Q.; Yang, G.; Kapila, V., Global output feedback tracking control of spacecraft formation flying with parametric uncertainty, Decision and Control, pp. 584 -589, Phoenix, AZ, December 7-10, 1999.

PAGE 183

165 [33]. Takashi Kida, Isao Yamaguchi, Yuichi Chida, and Takeshi Sekiguchi, On-orbit robust control experiment of flexible spacecraft ETS-VI, Journal of Guidance, Control, and Dynamics, Vol. 20 No. 5, pp. 865-873, September-October 1997 [34]. Richard J. Adams, and Siva S. Banda, Robust flight control design using dynamic inversion and structured singular value synthesis, IEEE Control Systems Technology, pp. 80-92, Vol. 1, June 1993. [35]. Datta N. Godbole, John Lygeros, and Shanker Sastry, Hierarchical hybrid control: an IVHS case study, IEEE Decision and Control Conference, pp. 1592 -1597, Vol. 2,Lake Buena Vista, FL, December, 1994. [36]. P.A.Stadter, Discrete event command and control for formation flying of distributed small spacecraft Systems, 13 th AIAA/USU Conference on Small Satellites, SSC-99-VI-4, Logan, Utah, August 23-26, 1999. [37]. M.Mesbahi, and F.Y.Hadaegh, Formation flying control of multiple spacecraft via graphs, matrix inequalities, and switching, Journal of Guidance, Control, and Dynamics, Vol.24, No.2, pp. 369-378, March-April 2001. [38]. Kang, W.; Sparks, A.; Banda, S, Multi-satellite formation and reconfiguration, American Control Conference, pp. 379-383, Vol. 1,Chicago, IL, June 28-30, 2000. [39]. Kang, W.; Xi, N.; Sparks, A, Theory and applications of formation control in a perceptive referenced frame, IEEE Decision and Control Conference, pp. 352-357 Vol. 1,Sydney, Australia, December, 2000. [40]. Andrew Robertson, Gokhan Inalhan, and Jonathan P. How, Formation control strategies for a separated spacecraft interferometer, American Control Conference, pp. 4142-4147, Vol. 6, San Diego, California, June 2-4, 1999. [41]. Johathan Lawton, Randal W.Beard, and Fred Y. Hadaegh, An adaptive control approach to satellite formation flying with relative distance constraints, American Control Conference, pp. 1545-1549, Vol. 3, San Diego, California, June 2-4, 1999. [42]. Ping Lu, Nonlinear Predictive Controllers for continuous systems, Journal of Guidance, Control, and Dynamics, Vol. 17 No. 3, pp. 553-560, May-June 1994 [43]. Haizhou Pan; Kapila, V., Adaptive nonlinear control for spacecraft formation flying with coupled translational and attitude dynamics, IEEE Decision and Control Conference, pp. 2057 -2062, Vol. 3, Orlando, FL, December, 2001. [44]. Gkhan Inalhan, Michael Tillerson, and Jonathan P. How, Relative dynamics and control of spacecraft formations in eccentric orbits, Journal of Guidance, Control, and Dynamics, Vol. 25, No. 1, pp. 43-60, January-February 2002.

PAGE 184

166 [45]. Gokhan Inalhan, Michael Tillerson, and Johathan P. How, Relative dynamics and control of spacecraft formations in eccentric orbits, Journal of Guidance, Control, and Dynamics, Vol.25, No.1, pp. 43-60, January-February 2002. [46]. Schiff, C.; Rohrbaugh, D.; Bristow, J. Formation flying in elliptical orbits, IEEE Aerospace Conference, pp. 37-47, Vol. 7, Big Sky, MT, March 18-25, 2000. [47]. Jasim Ahmed, Vincent T. Coppola, and Dennis S. Bernstein, Adaptive asymptotic tracking of spacecraft attitude motion with inertia matrix identification, Journal of Guidance, Control, and Dynamics, Vol. 21 No. 5, pp. 684-692, September 1998. [48]. Charng-Shi Wu, and Bor-Sen Chen, Adaptive attitude control of spacecraft: mixed approach, Journal of Guidance, Control and Dynamics, Vol. 24, No.4, pp. 755-767, July-August 2001. HH/2 [49]. Jongrae Kim, Jinho Kim, and John L. Crassidis, Disturbance accommodating sliding mode controller for spacecraft attitude maneuvers, AAS/GSFC 13 th International Symposium on Space Flight Mechanics, AAS 98-310, Greenbelt MD, May 11-15, 1998. [50]. Guang Q. Xing, and Shabbir A. Parvez, Nonlinear attitude state tracking control for spacecraft, Journal of Guidance, Control, and Dynamics, Vol. 24 No. 3, pp. 624-626, May 2001 [51]. James R. Wertz, Spacecraft attitude determination and control, Academic Publishers, Kluwer, 1978. [52]. John E. Prussing, Bruce A. Conway, Orbital Mechanics, Oxford University Press, Oxford, 1993 [53]. Roger R. Bate, Donald D. Mueller, and Jerry E. White, Fundamentals of Astrodynamics, Dover Publications, Inc., New York, 1971. [54]. Yunjun Xu, Norman Fitz-Coy, and Paul Mason, Coupled orientation and orbit maintenance in formation flight via sliding mode control, J. of British Interplanetary Society, Accepted [55]. Yunjun Xu, Norman Fitz-Coy, Generalized relative dynamics and control in formation flying system, 26th annual AAS Guidance and Control Conference, Breckenridge, Colorado, February 5 9, 2003. [56]. Alan Dow, Chairperson, Department of Mathematics, University of North Carolina at Charlotte, February 11, 2003, (last update), http://www.math.uncc.edu/~droyster/math3181 October 28, 2003, (last access).

PAGE 185

167 [57]. Richard H. Battin, An Introduction to the Mathematics and Methods of Astrodynamics, AIAA Education Series, New York, 1987 [58]. Andrew Robertson, Gokhan Inalhan, and Jonathan P. How, Spacecraft formation flying control design for the orion mission, AIAA Guidance, Navigation and Control Conference, pp. 1567-1576, Portland, Oregon, August, 1999. [59]. Miguel Rios-Bolivar, Adaptive canonical form for sliding mode control of uncertain nonlinear systems, 15 th Triennial World Congress of the International Federation of Automatic Control, Barcelona, Spain. [60]. Randal W.Beard, and Fred Y. Hadaegh, Fuel optmization for unconstrained rotation of spacecraft formations, The Journal of the Astronautical Sciences, Vol. 47, No. 3 and 4, pp.259-273, July-December 1999. [61]. Willem H Steyn, and Yoshi Hashida, An attitude control system for a low-cost Earth observation satellite with orbit maintenance capability, 13th AIAA/USU Conference on Small Satellites, SSC99-XI-4. [62]. David Folta, and David Quinn, A universal 3-D method for controlling the relative motion of multiple spacecraft in any orbit, NASA GSFC Flight Mechanics Symposium, pp. 405-406, Greenbelt, MD, May 1997. [63]. Gary E. Yale, and Brij N. Agrawal, Lyapunov controller for cooperative space manipulators, Journal of Guidance, Control, and Dynamics, Vo. 21, No. 3, pp. 477-485, May 1998. [64]. Seung-Bok Choi, and Myoung-Suk Kim, New discrete-time, fuzzy-sliding-mode control with application to smart structures, Journal of Guidance, Control, and Dynamics, Vol. 20, No. 5, pp. 857-865, September 1997. [65]. Mesbahi, M.; Hadaegh, F.Y., Mode and logic-based switching for the formation flying control of multiple spacecraft, American Control Conference, pp. 710-715, Vol. 2, Arlington, Virginia, June 25-27, 2001. [66]. J.F.Goulet, C.W.de Silva, V.J.Modi, and A.K. Misra, Hierarchical control of a space-based deployable manipulator using fuzzy logic, Journal of Guidance, Control and Dynamics, Vol. 24, No.4, pp. 395-405, July-August 2001. [67]. Daniele Mortari, Euler-q algorithm for attitude determination from vector observations, Journal of Guidance, Control, and Dynamics, Vol. 21, No. 2, pp. 328-335, Issue Mar 1998. [68]. Rick Lind, EML 5311 Control Theory II, Class notes, Spring 2003, University of Florida [69]. Gary J. Balas, John C. Doyle, Keith Glover, Andy Parkard, and Roy Smith, Mu Analysis and Synthesis Toolbox, Mathworks Inc., Natick, MA. 2001.

PAGE 186

168 [70]. J.E. Slotine and W. Li, Applied Nonlinear Control, Prentice Hall, Englewood Cliffs, New Jersey. [71]. Sergey Drakunov, and Vadim Utkin, Sliding mode observers. Tutorial, 34 th Conference on Decision and Control, pp. 3376-3378, Vol. 4, New Orleans, LA December 1995. [72]. Yunjun Xu and Norman Fitz-Coy, Formation flight system via robust control, 27th annual AAS Guidance and Control Conference, Breckenridge, Colorado, February 4 8, 2004. [73]. Sergey Edward Lyshevski, Sliding modes and soft switching control in dynamics systems, American Control Conference, pp. 646-650, Vol. 1, Chicago, Illinois, June 2000. [74]. G. Monsees, J.M.A. Scherpen, Adaptive switching gain for a discrete-time sliding mode controller, American Control Conference, pp. 1639-1643, Chicago, Illinois, June, 2002. [75]. B. W. Bekit, J. F. Whidborne, and L. D. Seneviratne, Sliding mode control for robot manipulators using time-varying switching gain and boundary layer, UKACC International Conference on Control, pp34-40, Swansea UK, September 1-4, 1998. [76]. K. David Young, Vadim I. Utkin, Umit Ozguner, A control engineers guide to sliding mode control, 1996 IEEE Workshop on Variable Structure Systems, pp. 1-14, Tokyo, Japan, December 5-6, 1996. [77]. V. I. Utkin, Sliding mode control in dynamic systems, Proceedings of the 32 nd Conference on Decision and Control, pp. 2446-2451, Vol. 3, San Antonio, Texas, December 1993. [78]. Vadim Utkin, Sliding mode control in mechanical systems, Industrial Electronics, Control and Instrumentation, pp. 1429-1431, Vol. 3, Bologna, Italy, September 5-9, 1994. [79]. Yunjun Xu, Norman Fitz-Coy, and Paul Mason, Thruster limitation considerations for formation flight control, Mission Engineering and Systems analysis Division, Flight Mechanics Symposium, GSFC, Greenbelt, MD, October 28-30, 2003. [80]. Jing Wang, and Yuanyuan Zhao, The optimization of PI controller parameters using genetic algorithm in the DC speed control system, 3 rd World Congress on Intelligent Control and Automation, pp. 545-548, Hefei, P.R. China, June 28 July 2, 2000.

PAGE 187

169 [81]. Ta-Tau Chen, and Tzuu-Hseng S. Li, Integrated fuzzy GA-based simplex sliding-mode control, International J ournal of Fuzzy Systems, Vol. 2, No. 4, December, 2000. [82]. Powell, M.J.D., "A fast algorithm for nonlinearly constrained optimization calculations," Numerical Analysis ed. G.A. Watson, Lecture Notes in Mathematics Vol. 630, Springer Verlag, NY, 1978. [83]. Yunjun Xu, Norman Fitz-Coy, and Paul Mason, Genetic algorithm based sliding mode control in the leader/follower satellites pair maintenance, 2003 AAS/AIAA Astrodynamics Specialist Conference, Big Sky, Montana, August 3-7, 2003. [84]. Yunjun Xu, Norman Fitz-Coy, and Paul Mason, Maintaining and reconfiguration control strategy in formation flight system , J. of British Interplanetary Society, Accepted. [85]. Wei Kang, Andy Sparks, and Siva Banda, Coordinated control of multisatellite systems, Journal of Guidance, Contro l, and Dynamics, Vol. 24, No. 2, pp. 360-369, March-April 2001.

PAGE 188

BIOGRAPHICAL SKETCH Yunjun Xu was born in Jiangsu, China, P.R., on October 30, 1973. In July 1996 and April 1999, he received the Bachelor and Master of Science from Nanjing University of Aeronautics and Astronautics, Jiangsu, China, P.R., respectively. In August 1999, he became a graduate student in the Department of Mechanical and Aerospace Engineering at the University of Florida with a graduate teaching and research assistantship. On January 7, 2002, a lovely son, Ronald Xu, was born to Yunjun Xu and his wife Yanli Xia. On May 2002, he received the Master of Science from the Department of Electrical and Computer Engineering at the University of Florida as a concurrent degree. He is currently completing the requirements for the degree of Doctor of Philosophy in the Department of Mechanical and Aerospace Engineering. His research interests are linear and nonlinear control theory and application, orbit and spacecraft attitude dynamics and control, and multi-agents system. 170


Permanent Link: http://ufdc.ufl.edu/UFE0002318/00001

Material Information

Title: Dynamics and Control for Formation Flying System
Physical Description: Mixed Material
Copyright Date: 2008

Record Information

Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved by the source institution and holding location.
System ID: UFE0002318:00001

Permanent Link: http://ufdc.ufl.edu/UFE0002318/00001

Material Information

Title: Dynamics and Control for Formation Flying System
Physical Description: Mixed Material
Copyright Date: 2008

Record Information

Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved by the source institution and holding location.
System ID: UFE0002318:00001


This item has the following downloads:


Full Text











DYNAMICS AND CONTROL FOR FORMATION FLYING SYSTEM


By

YUNJUN XU


















A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL
OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT
OF THE REQUIREMENTS FOR THE DEGREE OF
DOCTOR OF PHILOSOPHY

UNIVERSITY OF FLORIDA


2003































Copyright 2003



by
Yunjun Xu





























Dedicated to

my parents,
who always loved;
my wife, Yanli
who always encouraged;
my son, Ronald
who always delighted.















ACKNOWLEDGMENTS

I would like to express my gratitude to my advisor, Dr. Norman Fitz-Coy, for his

support, patience, and encouragement throughout my graduate studies. It is not often that

one finds an advisor and colleague that always find the time to listen to the little problems

and roadblocks that unavoidably crop up in the course of performing research. His

technical and editorial advice was essential to the completion of this dissertation and has

taught me innumerable lessons and insights on the workings of academic research in

general.

My thanks also go to the members of my committee, Drs. Andrew Kurdila, Rick

Lind, Peter Ifju, and Michael Nechyba, for reading previous drafts of this dissertation and

providing many valuable comments that improved the presentation and contents of this

dissertation. I would like to thank Dr. Paul Mason for his help during my whole research

period.

The friendship of colleagues in the Dynamics and Control Group is much

appreciated and has led to many interesting and good-spirited discussions relating to this

research.

Last, but not least, I would like to thank my wife, Yanli, and my son, Ronald, for

their understanding and love during the past few years. Their support and encouragement

were in the end what made this dissertation possible. My parents receive my deepest

gratitude and love for their dedication and the many years of support that provided the

foundation for this work.















TABLE OF CONTENTS

page

ACKNOW LEDGM ENTS .................................................................... iv

LIST OF TABLES................... .................. .................. .. ......... viii

LIST OF FIGURES................... .................... .... ............................ix

ABSTRACT ................. ............. ............................... ..xvii

CHAPTER

1 INTRODUCTION ................ ...... ....................... ............... 1

1.1 M otivation............................... ........... ......... 1
1.2 R research Scope................................. ............... .... ........ 8
1.3 Dissertation Outline............................................ ......... ..9

2 MULTIPLE VEHICLES DYNAMIC MODEL.......................................11

2.1 Introduction ............. ....... ........... ..... ............ ...... .... .. 11
2.2 Discussion Of Multiple Vehicles Dynamic Model...................... 14
2.3 General 2-Body Relative Motion Models....................................17
2.3.1 Relative Motion Dynamic Model (Point Mass Case)............18
2.3.2 Simulation Results and Discussion (Point Mass Case)..............28
2.3.2.1 Initial condition ................. .................. .................28
2.3.2.2 Results and discussion.........................................28
2.3.3 Specialized Relative Motion Model (Point Mass Case)...........34
2.3.3.1 Linear model (Clohessy-Wiltshire equation) ..............34
2.3.3.2 Circular reference orbit nonlinear model.................. 36
2.3.3.3 Elliptic reference orbit nonlinear model....................37
2.3.4 Simulation Results and Discussion (Point Mass Case)...............38
2.3.5 Relative Motion Dynamic Model (Rigid Body Case)................40
2.3.5.1 D esired attitude ............ ... ............................. 40
2.3.5.2 Attitude m odel............ .. .......... ... ...........47
2.3.6 Simulation Results and Discussion (Rigid Body Case)..............49
2.4 P erturbation ....................................................... ....... 51
2.4.1 A tm ospheric D rag................................. .......... .......52
2.4.2 Earth Oblateness......................................................... 53















2.4.3 Perturbation Effects on the Generalized Dynamics Model........56
2.4.4 Perturbation Effects on Keplerian Elements.........................58
2.5 Conclusion...................................64

3 CONTROL STRATEGIES: ORBIT............. .................................. 65

3.1 Introduction ................ ............... ................. .......... 65
3.2 Performance Criterion Of The Leader-Follower Pair......................69
3.3 Relative Position Controller Design ................ ... ........................70
3.3.1 Linear-Quadratic State Feedback Regulator (LQR)...............70
3.3.1.1 LQR ........... .......................................... 70
3.3.1.2 LQR synthesis m odel..........................................71
3.3.1.3 L Q R analysis m odel.............................................74
3.3.2 H //p Controller ......... ..... ................ .... ...............76
3.3.2.1 Structured singular value (/u) theory .......................76
3.3.2.2 / synthesis m odel................................... .... 77
3.3.2.3 /u analysis nodel............................. ....... .......... 81
3.3.2.4 N P, R S, and RP ............. ................. ............ 81
3.3.3 Feedback Linearization ............... .........................83
3.3.3.1 Feedback linearization........................... .................83
3.3.3.2 Feedback linearization analysis model ........................89
3.3.4 Sliding Mode Controller (SMC).....................................90
3.4 Relative Position Simulation Results and Discussion.......................92
3.4.1 State Feedback Gain ............... ..................... ............93
3.4.2 Robustness and Performance Analysis............ ............ 95
3.4.3 Time Response Based on C-W Model................... ............ 96
3.4.4 Time Response Based on C-W Model with Uncertainties...........98
3.4.5 Time Response Based on Actual Nonlinear Model with Thruster
Lim itation ............. ...................... ................... ....... 101
3.4.6 Conclusion...................... ................ 102

4 CONTROL STRATEGIES: ATTITUDE AND ORBIT ...........................104

4.1 Introduction ........................................................ ........... 104
4.2 Relative Attitude Control................... ................ .................104
4.2.1 Control Architecture ............ ......... .......... 104

4.2.2 The Inner Loop Controller ......... .............. ................106
4.3 Combined Relative Position and Attitude Controller .....................112
4.4 Boundary Layer Control................... ............... ... .............. 116















4.4.1 Hard Switching Control........... ........... ..............117
4.4.2 Boundary Layer Control............... .......... ............... 117
4.4.3 Sim ulation R esults................... ................ ....... ........ 118
4.5 G enetic A lgorithm ..................................... ........... ........ ...125
4.5.1 Optimization Problem ............. ................. .. ........ 125
4.5.2 Genetic Algorithm................................. ...............127
4.5.3 Sim ulation R results ................... ................. .... ..........131

5 RECONFIGURATION CONTROLLER IN FORMATION FLYING........... 136

5.1 Introduction ............. ........... .... ............................ 136
5.2 Formation Reconfiguration Strategy .......................................140
5.3 Perceptive Frame Control (PFC)................ .............................143
5.4 Controller D esign................... ......................... ................145
5.4.1 Perceptive Frame Controller ...... ............................ 145
5.4.2 Formation Maintenance Controller...............................146
5.4.3 Formation Reconfiguration Controller........... ...............147
5.4.4 SIMULINK Block Diagram............. ......... ......... 147
5.5 R results ............. ................... ...................... ... .... 149
5.6 Sum m ary...................................................... ......... 156

6 CONCLUSION AND FUTURE WORK.............................. ............157

6 .1 C onclu sion ................................................ 157
6.2 Future Work............................................... .... .........158

APPENDIX

QUATERNION AND LIE DERIVATIVE................................................160

A 1 Quaternion.................................................. ..... ...... 160
A.2 Lie Derivative................... ...........................................161

LIST OF REFERENCES.................................................................. 162

BIOGRAPHICAL SKETCH........... ............................................... 170















LIST OF TABLES

Table page

2-1 Different cases been simulated............. ......... .................... 29

2-2 Classical elements for an Earth satellite........... ............ ................59

3-1 G ain and phase m argins............... ................ .............. ....... ......95

3-2 Robustness performance and controller order comparison.......................96

3-3 Fuel consum ption................................... ..................... .......... 98

3-4 Fuel consum ption............................... ................ ..............101

4 -1 C ases ................... ................... ..... .................. .... ....... 107

4-2 C ases ................... ................... ..... ................... .... ....... .110

4-3 C ases ............. ............................. .................... ... ......... 113

4-4 GA parameters' definition and the data used..................................... 132

4-5 Simulation results of GA + SMC + BLC............... .......... ..............133

5-1 The initial orbit information of the formation satellites..................... 139

5-2 Orbit information of R1 and R2 in figures 3 and 4 ........................... 142

5-3 Formation maintenance controllers for follower 1 and follower 2.............147

5-4 Formation reconfiguration controllers for redundant 1 and redundant 2........147

5-5 Configuration code for in-plane formation flight simulation........ ......... 149















LIST OF FIGURES

Figure page

1.1 Small satellite for DS1 ............................................ .... ...........

1.2 The probe and m mission outline ................. .............. .............................. ............... 5

1.3 T h e satellite for E O -1 ......... .. ....................................................... ........ .......... .. .. .6

1.4 F orm action flying of ST 5 ................................................................ ........................ 7

1.5 Size of the micro-satellite and the communication architecture ...............................7

2.1 Classification of the relative dynamics model ...................... .................... ........... 11

2.2 N bodies system ..................................... .......................... ..... ........... 15

2 .3 Z ones of variations........ ........................................................................... ....... ... ... 16

2 .4 R elative dynam ic m odel ........................................................................ ...................18

2.5 Coordinate system of the relative position model......................................................19

2.6 Coordinate system of the leader satellite's frame....................................................19

2.7 Orbital elements ..................................................................... ........ 25

2.8 O rbital inform action ......................................................... ........... ............ .... 1

2.9 R elative position in x direction ....................................................................... 31

2.10 R elative position in y direction .......................................................... ............... 31

2.11 R elative position in z direction ............................ ..... ......... .......................... 31

2 .12 R elativ e position in 3D ........................................................................ ...................32

2.13 M odel difference ........ ....... ...... ................................. ..... 32

2 .14 C ase 9 m odel difference...................................................................... ...................32

2 .15 C ase 20 m odel difference...................................................................... ..................34



ix
















2.16 R elative position in x direction ........................................ ............... ............... 39

2.17 R elative position in y direction........................ ... .......................... ............... 39

2.18 Difference between CW and generalized model .................................. ...............39

2.19 Satellite attitude body fram e .............................................. ...................... ............... 41

2.20 Spherical right triangle.............................................. .................... ............... 42

2.21 Desired attitude I..................... ........................... .......43

2.22 D desired attitude II .......................... .......... .. .......... ...... ....... 43

2.23 D desired attitude III .................. .............................. ...... ................. 44

2.24 Error angles in the coordinate system. EL is the elevator error angle, whereas, AZ
is the azim uth error angle ...................... .. .. ............................... .... ........... 46

2.25 Theta 1 (EL) and Theta2 (AZ) in radians ........................................ ............... 50

2.26 Theta 1 (EL) and Theta2 (AZ) in radians ........................... ...................... 50

2.27 Thetal (EL) and Theta2 (A Z) ....................................................... ... ............... 51

2 .28 A tm sphere drag .......................... ...................................................... 53

2.29 The perturbation stems from the Earth oblateness.............................. ..............54

2.30 A tm sphere drag effects ....................................................................... ..................57

2 .3 1 J2 perturb action effects ......................................................................... ....................57

2.32 Combined perturbation effects................. ........... .... ............... 57

2.33 Error in position ................ ........................... ................ 58

2.34 Error in attitude ....................... ................................... ......... ......... 58

3.1 C om bined state space............ ............................................................... .......... ....... 72

3.2 Analysis model of LQR controller without uncertainty ...........................................75

3.3 Analysis model of LQR controller with uncertainty...... ..........................................75


x















3.4 Typical loop gain and phase m argin m odel ........................................ .....................76

3.5 Open loop model with parametric uncertainty..................................... ............... 79

3.6 Synthesis M odel ................................ .. ................................ .......... 80

3.7 Synthesis m odel input-output relation ........................................ ........................ 81

3.8 Analysis m odel for nom inal system ....................................................... 81

3.9 Analysis m odel for actual system ............................................................................81

3.10 Closed-loop model with uncertainty..................... .. ................ ................. 82

3.11 Equivalent model 1 ................................. .. .. ........................ 82

3.12 Equivalent m odel 2 ................................... ............. .. ............82

3.13 Closed-loop analysis model based on the feedback linearization.................................89

3.14 Closed-loop analysis model for the linear system ................................ ...............90

3 .15 R elativ e position error in x ........................................ .............................................96

3.16 R elative position error in x ...................................................... ... ............. 96

3.17 R elative position error in y .................................................. ......................... 97

3.18 R elative position error in y ......................................................................... ............ 97

3.19 R elative position error in z............................ .............. .................. ............... 97

3.20 R elative position error in z................................................. ............................... 97

3.21 Control force in x ........................................ ..................... .... .... ............... 98

3.22 C control force in y .................. .................................................. .... ...... 98

3.23 C control force in z .......................... ...... ................ .. .. ............. ........ 98

3.24 R elative position error in x ........................................ .............................................99

3.25 R elative position error in x ........................................ .............................................99















3.26 R elative position error in y ........................................ .............................................99

3.27 R elative position error in y ........................................ .............................................99

3.28 R elative position error in z................................................. ............................... 100

3.29 R elative position error in z................................................. ............................... 100

3.30 Control force in x ................................ ......................... .. ............ 100

3.31 C control force in y ................ ............................ .................... ....... 100

3.32 Control force in z ................................. .. .. ................... ....... 101

3 .33 R elativ e position error in x ............................................ ......................................... 10 1

3 .34 R elativ e position error in y ............................................ ......................................... 10 1

3.35 R elative position error in z................................................. ............................... 102

4 .1 A attitude control architecture ........................................ .............................................105

4.2 Relative position difference between exact and uncertainties model...........................108

4.3 Relative position difference between exact and uncertainties model...........................108

4.4 Elevator error angle difference between the exact and uncertainties models...............09

4.5 Azimuth error angle difference between the exact and uncertainties models ................109

4.6 Elevator error angle if no control is applied ..................... ....... ............. ............... 109

4.7 Azimuth error angle if no control is applied................. ..................109

4.8 Elevator error angle under the sliding mode controller .............................................111

4.9 Azimuth error angle under the sliding mode controller in transient stage from 0 600
seconds ................................................................ .......... .......... 111

4.10 Azimuth error angle under the sliding mode controller from 5000 5800 seconds ....111

4.11 E elevator error angle in case 5............................................... ............................. 111

4.12 Azim uth error angle in case 5 ...................................................................... 112















4.13 R elative position in x direction .............................................. ................................. 113

4.14 R elative position in y direction ....................................................................... 113

4.15 Relative position in z direction ..................................................... ................... 114

4 .16 T he elevation error angle ................................................................ .......................114

4.17 Zoomed in elevation error angle ..................................................... .................. 114

4.18 Zoomed in azimuth error angle ....................................................... ...............114

4.19 Zoomed in azimuth error angle ....................................................... ...............114

4.20 Force com m ands in x direction........................................................................... ... 114

4.21 Zoomed in force commands in x direction ......................................... ...............115

4.22 Force com m ands in y direction.............................................. ........................... 115

4.23 Zoomed in force commands in z direction ........... ..............................................115

4.24 T orque com m hands in axis 1 ............................................................................ ..... 115

4.25 Torque com m hands in axis 2 ........................................................................... ...... 116

4 .26 Sliding surface in SM C ............................................... .................... .......................118

4.27 Sliding surface in SMC + BLC ....................................................... ...............118

4.28 Relative position in x direction, transient stage .....................................................119

4.29 Relative position error in x direction, stable stage.......................................................119

4.30 Relative position in y direction, transient stage .....................................................119

4.31 Relative position error in y direction, stable stage.......................................................119

4.32 Relative position in z direction, transient stage. .................. .............. ............... 119

4.33 Relative position in z direction, stable stage ............. ........................119

4.34 Sliding surface for the relative position in x direction, transient stage.........................120















4.35 Sliding surface for the relative position in x direction, stable stage ...........................120

4.36 Sliding surface for the relative position in y direction, transient stage.........................120

4.37 Sliding surface for the relative position in y direction, stable stage ...........................120

4.38 Sliding surface for the relative position in z direction, transient stage.........................120

4.39 Sliding surface for the relative position in z direction, stable stage. ..........................120

4.40 The elevator error angle in the transient stage............. .............................................. 122

4.41 The elevator error angle in the stable stage ........................................ .....................122

4.42 The azimuth error angle in the transient stage..........................................................122

4.43 The azimuth error angle in the stable stage ...................................... ............... 122

4.44 Force commands in x direction of the local frame, in transient stage .........................123

4.45 Force commands in x direction of the local frame, in stable stage............................. 123

4.46 Force commands in y direction of the local frame, in transient stage .........................123

4.47 Force commands in y direction of the local frame, in stable stage............................. 123

4.48 Force commands in z direction of the local frame, in transient stage..........................124

4.49 Force commands in z direction of the local frame, in stable stage............................. 124

4.50 Torque commands in body axis 1, in transient stage...................................................124

4.51 Torque commands in body axis 1, in stable stage ................................ ............... 124

4.52 Torque commands in body axis 2, in transient stage...................................................124

4.53 Torque commands in body axis 2, in stable stage ................................ ............... 124

4.54 A generation of the chrom osom e...................................................................... ...... 128

4.55 Evaluation, fitness, and selection ................................................... ...............129

4.56 Crossover .................................... ................... ................. ......... 130















4 .5 7 M u tatio n .............................. ................................................................. ............... 13 0

4.58 The relative position error in x-direction. ................................. ............................... 134

4.59 Control force in x-direction. ............................................... .............................. 134

4.60 The relative position error in y-direction. ................................. ............................... 134

4.61 Control force in y-direction. ............................................... .............................. 134

4.62 The relative position error in z-direction. ........................................ ............... 134

4.63 C control force in z-direction ......... ................. ................... .................. ............... 134

4.64 Elevation error angle ........ ................................. ........ ............... .....135

4.65 Torque in body axis 1. .. ................ .......................................... ........ .. ............... 135

4.66 Azimuth error angle. ...................................... ........................135

4.67 Torque in body axis 2. .................................................... ............ .. 135

5 .1 A b solute ......... ..... ........... ................ ........ ................... ............... 137

5.2 Leader-follow er Structure................................................................. ............... 137

5.3 A b solute and interconnection ........................................... ........................................ 137

5 .4 Su b form action ............................. ........................................................... ............... 13 7

5.5 3-satellite formation architecture with two redundant satellites ................................139

5 .6 If F 2 fails ......... ..... .......... .............. .....................................................14 1

5 .7 If F 1 fails .............................................................................. 14 2

5.8 Perceptive frame controller's structure ................................. ...............145

5.9 SIM U LINK program diagram ............................................... ............................. 148

5.10 Relative position in x local frame. (a) whole range, and (b) from 5670 7500
second ........................................................................... 150




xv















5.11 Relative position in y local frame, (a) whole range, and (b) from 5500-10000
se c o n d s ...................................... ............................... ................ 1 5 0

5.12 Relative position differences in z local frame, (a) 0-800 seconds, and (b) 5650-6500
seconds ................................................................ ... ..... ......... 151

5.13 Force commands in x local frame. (a) whole range, and (b) 5600-10000 seconds ......152

5.14 Force commands in y local frame. (a) whole range, and (b) 5600-9000 seconds ........152

5.15 Force commands in z local frame. (a) 0-800 seconds, and (b) 5900-6300 seconds .....153

5.16 Elevation error angle. (a) whole range, and (b) 5670-15000 seconds ..........................154

5.17 Azimuth error angle. (a) whole range, and (b) 5900-7000 seconds............................154















Abstract of Dissertation Presented to the Graduate School
of the University of Florida in Partial Fulfillment of the
Requirements for the Degree of Doctor of Philosophy

DYNAMICS AND CONTROL FOR FORMATION FLYING SYSTEM
By

Yunjun Xu

December 2003


Chair: Norman Fitz-Coy
Major Department: Mechanical and Aerospace Engineering

Due to the potential advantages of the micro-satellites formation flying system,

such as cost reduction, flexibility, observational baseline increase, and reliability

improvement, numerous formation related missions have been proposed and seriously

studied. Along these many potential benefits, come a variety of challenges. Dynamics

and control of the multi-satellites in the formation is one of these issues and is addressed

in this dissertation. To date, most formation flying dynamics and control studies are

based on the linear Clohessey-Wiltshire equations and none address the problem with a

coupling between the attitude and orbit motions. Furthermore, commercial software, such

as STK and FreeFlyer, is not capable of feedback control simulation. These observations

motivate this research to develop a generalized coupled attitude-orbit relative motion and

to design controllers based on this model.

The developed generalized dynamics model is based on the leader-follower

structure and the desired attitude is coupled with the real time relative position due to the

formation mission requirement. This nonlinear model is applicable to any closed


xvii









reference orbit cases and incorporates a coupling between the attitude and orbit. Also,

atmospheric drag and J2 perturbations are included.

The formation flying control system can be regarded as a two-tiered architecture.

In the lower level, continuous controllers are designed to meet the objectives of each

vehicle, and in the higher level, discrete time event determines the formation

configuration and reconfiguration. After comparison with LQR, H, / / synthesis, and

feedback linearization strategies, sliding mode control is chosen as the basic control

method for the lower level in this dissertation. Furthermore, boundary layer control and

genetic algorithm are implemented in order to reduce chatters with the implementation of

the sliding mode control and reduce the fuel consumption. In the upper level, perceptive

frame concept is adopted for the discrete event determination. It is demonstrated that the

proposed control architecture is capable of tracking and maintaining the desired position

and attitude under different formation missions.


xviii














CHAPTER 1
INTRODUCTION



1.1 Motivation

In the 1950's, the first satellite was launched by the Soviet Union [1]. From that

moment on, more and more satellites have been launched into space for the purpose of

Earth observation, deep space exploring, military surveillance, navigation (GPS),

commercial and military communication, and weather prediction.

To date, conventional monolithic satellites typically have the following

deficiencies. The first deficiency is their large size and/or weight. A significant portion of

a mission's budget is spent on a launch vehicle. The larger a satellite, the larger and more

costly a launch vehicle is required to boost it into space. The second deficiency of the

conventional satellite is its lack of flexibility for missions. When the satellite is located in

its orbit, the mission is fixed. The instruments of the satellite will limit the observing

baseline. The last deficiency is the bad redundancy. Whenever a failure occurs on the

satellite, typically the entire mission fails or the redundant subsystem may need to be

used, which will increase the payload of the satellite.

Therefore, a new concept, formation flying system or distributed space system,

was introduced in the early 80's by Hartman and Hadaegh et al. [2, 3]. In 1998, the

formation flying system has received renewed interest by numerous government

agencies. The formation flying systems are multiple mini, micro, or nano satellites that

can work together in order to perform a unified mission [4]. Such a system will take the









place of the conventional monolithic satellites. "This system can range from global

constellations offering extended service coverage to clusters of highly coordinated

vehicles that perform distributed sensing." [4].

NASA's New Millennium Program [5] proposes the properties of the formation

flying system as the following:

+ Fleets of vehicles

+ Vehicles interact and cooperate with one another to achieve mission goals

+ Vehicle pointing managed collectively

+ Vehicle positioning managed collectively

+ Fleets evolve over time to extend and expand capability

+ Self-controlling vehicles eliminate the need for extensive ground support

With respect to the conventional single large satellites, a formation flying system

may provide numerous advantages [5-8]. Some of these advantages are summarized

below.

+ Cost Reduction: Due to the size and weight reduction, the launch cost will decrease
substantially. Using the Space Shuttle as an example, the cost to deploy a payload into a
low-Earth orbit is about $13,200 per kilogram ($6,000 per pound) [5]. Most current
satellites have masses on the order of 1,000 Kg. But in formation flying system, the mass
of a micro-sat will be about 100 Kg, and that of a nano-sat will be approximately 20 Kg.
Therefore, at least 50 nano-satellites can be launched into the space at the same cost with
potentially more powerful functions.

+ More flexible, efficient use of resources: The formation of the multiple satellites can be
changed autonomously or manually according to different missions or error conditions.

+ Extensive co-observing programs can be conducted without using extensive ground
support. For instance, in the leader-follower communication architecture, the satellites do
not need to communicate with the ground station all the time. Only the leader satellite is
in charge of the task and the autonomous handling software will determine when the
communication with the ground support is necessary.









+ Increase precision and observational baselines. For example, "A team of powerful X-
Ray telescopes unlike anything this part of the universe has ever seen. These telescopes
will work in unison to simultaneously observe the same distant objects, combining their
data and becoming 100 times more powerful than any single X-ray telescope that has
come before it."[5]

+ Enhanced survivability and increased reliability: In a case in which a certain number of
the satellites in the formation system failed, the mission may still be accomplished. This
is an important advantage over the conventional monolithic satellite.

With these potential benefits, however, come a variety of technical challenges.

Some of them are the following:

1. Almost all space missions are flown with a single spacecraft. Controlling spacecraft or
satellites in formation flight is very complicated.

2. Hybrid system and reconfiguration control is required.

3. The ground-to-multiple-satellite communications path is more complex than the
ground-to-single-satellite communications.

4. Because of the limitation in size and weight, designing and manufacturing the mini,
micro, or nano size components for the satellites in formation flying system are not easy
issues [9].

Although complete formation flying systems have not been implemented in space,

NASA and other government agencies have begun to launch several missions in order to

validate critical technical issues, and to increase the confidence of this brand new space

exploration era. NASA's New Millennium Program, which includes such missions as

Deep Space 1, Deep Space 2, Earth Observing 1, Earth Observing 3, and Space

Technology 5 [5, 10-12], is an example of a government program that is specific to

the spacecraft formation. These missions will be briefly discussed.

The Deep Space 1 (DS1) mission was launched on October 24, 1998. The retired

date of DS1 was December 18, 2001. Figure 1.1 shows the picture of the small satellite.

The total mass of DS1 is 486.3 Kg. Several techniques for the multiple satellites have









been tested in this mission; for example, low-mass communications devices, which are

half as heavy and 1/3 as costly as before, microelectronics and spacecraft structure, which

tried to demonstrate futuristic technologies for making the spacecraft smaller, lighter, and

more efficient, and autonomous operations system that can plan, make decisions, and

operate by itself.












Figure 1. 1 Small satellite for DS1



The Deep Space 2 (DS2) project includes two probes, piggybacked aboard the

Mars Polar Lander. It was launched on January 3, 1998. The two probes were slammed

into Martian soil on December 3, 1999. Figure 1.2 shows the probe and the mission

diagram. DS2 is aimed to design, develop, and flight validate advanced systems and

instruments in deep space in order to lower the cost and risk for future formation

missions. The total mass of the probe is 6.5 Kg [13].
























Figure 1. 2 The probe and mission outline



Although it is in still in an elementary stage, EO-1 is actually the first prototype in

the New Millennium Program of a formation flying system. EO-1 works with the

LandSat-7 to form a simple formation once or twice each day. Landsat 7 is a U.S.

satellite used to acquire remotely sensed images of the Earth's land surface and coastal

regions. At this point, EO-1 and LandSat-7 will image the same point in order to get the

piecewise continuous data. The EO-1 small satellite's orbit is a 705 Km circular, sun-

synchronous orbit inclined at a 98.70 to the equator. It is still in an elementary stage of the

formation-flying concept because EO-1 can only coordinate satellites when they are in

one or two points of the orbit and the complicated technologies, such as multiple

vehicles' communication, fault tolerant algorithm, and autonomous system, have not

been tested. The mission started on November 21, 2000. Figure 1.3 shows a picture of the

EO-l's small satellite. The mass of the small satellite is 90 Kg.


























Figure 1. 3 The satellite for EO-1



The objective of the EO-3, which is also called GIFT-IOMI [14], is to improve

the observation of the temperature, moisture, wind velocity, and test the new remote

sensor, which will be used in future missions. The planned launch date is 2005-2006.

The mass of the instrument is about 100 Kg.

Satellite Technology 5 (ST-5) is composed of three miniature satellites with a

mass of 19.4 Kg and 54.2 cm X 28.6 cm dimensions. Although it is small, it has full

"service" including guidance, navigation, control, propulsion, bandwidth, and

communication systems. This mission will be launched in 2004.
























Figure 1. 4 Formation flying of ST5
iLmL

























I ....................................
54.2 cm(21.3")
Figure 1. 5 Size of the micro-satellite and the communication architecture.



The right-hand side of Figure 1.5 shows the communication architecture that will

be used in this mission. Every individual satellite will communicate with the ground base.

In addition to the NASA related projects, there are others working in this area, for

example, Constellation X-ray, ACE [15], UoSAT-12, LISA and MAXIM. Furthermore,

there are some other projects related to the formation flight systems in other application,

such as TechSat 21 [16] of DoD and Air Force formation flight system [17].









A formation flying system is a complicated comprehensive system with numerous

sub-systems. The ORION micro-satellite is a correct example of an element of a

formation flying system. Its power, communication, dynamics and control chip, GPS

receiver, structure, orbit route, orbit condition, shape of the formation, and sensor are

discussed in Ref. [18]. In this dissertation, it is impossible to cover all these topics. In the

next section, the specific research topics of a formation flying system addressed in this

dissertation will be expressed.



1.2 Research Scope

In this section, the specific research topics that are covered in this dissertation will

be summarized, and a brief introduction of each is presented. A detailed literature review

of each topic is presented in the individual chapter.

Although the component technologies may be ready, the management and control

of a distributed constellation is still a major area of concern. Therefore, the dynamic

models and the controllers are two of the key issues discussed in this dissertation. The

literature abounds with papers that use the Clohessey-Wiltshire (C-W) equations as the

dynamic model upon which the controller is designed [19-30]. C-W equation, however, is

linear relative motion model that is restricted to the near circular orbits, and thus neglects

a large class of orbit cases. Recently, some papers [30, 31] studied the nonlinear relative

motion problem, but are still restricted to the circular reference orbit case. In this

dissertation, the dynamic models are extended to the generalized, non-circular, closed

orbits and are not linearized. These models include gravitational perturbations due to

oblateness of the central body as well as atmospheric drag.









The design of the controllers for the leader-follower hierarchical structure of

satellite formation is also a key technology which needs to be solved. To date, most

designs are based on the PID or LQR controller methodology. Recently, some elements

of robust control have been utilized in the controller design [4, 32-34]. In this thesis,

nonlinear controllers based on sliding mode control theory will be utilized. Various

controller designs will be evaluated based on performance criteria that address the fuel

consumption (proportional to magnitude of control effort) and the errors between the

actual and desired states.

Due to planed mission or unplanned changes (e.g., satellite failure), the formation

may need to be reconfigured. A hybrid control architecture will be used to address the

reconfiguration issue. Hybrid control methodology includes a multi-tiered control

architecture that has been used in the smart highway system [35]. In Stadter [36], the

high-level tier uses a discrete even command to achieve formation reconfiguration.

Graphs, matrix inequalities, and switching are also methods to deal with the system level

problem [37]. In this dissertation, a perceptive frame concept will be utilized and it will

be shown to have a good and smooth performance [38, 39].

1.3 Dissertation Outline

This dissertation is organized as follows. In Chapter 1, the basic concept and the

advantages of the satellite formation flying system were addressed. The state-of-the-art of

the formation-flying projects was illustrated. Also the research scope of this dissertation

will be provided.

Chapter 2 addresses the detailed dynamic models needed to analyze a spacecraft

formation flying system. The existing models in the literature will be presented and

discussed; their deficiencies will be high lighted. After that, the generalized nonlinear









dynamic model of the leader-follower satellites pair will be derived from the first

principles. It will be shown that the specialized models that appeared in the literature can

be derived from the generalized model under proper conditions. Also, in this chapter, the

combined attitude and orbit model will be developed. Finally in Chapter 2, the

perturbation terms and their effects on the classical six-orbital elements will be discussed.

Chapters 3, 4, and 5 address the control strategies for the formation flying system.

A two-tiered control structure is utilized in the formation system. The "lower" tier

addresses the specific vehicle control actions and the "upper" tier addresses the hierarchy

and reconfiguration control actions. The Chapters 3 and 4 will be focused on designing

the lower part, which is the controller for the relative position and attitude's maintaining

and tracking for the leader and follower satellites pair. The controller study will be based

on sliding mode control methodology, and the performance will be compared with the

LQR, p synthesis, and feedback linarization methods, which often appeared in the

literature. All controllers designed are rigid spacecraft approach. Furthermore, some

improvements for the sliding mode control are discussed. Chapter 5 addresses the higher

tier of the controller and focused on reconfiguration of the whole formation flying

system. The controller design is based on the perceptive frame concept. It will be shown

that the overall system is stabilized and controllable.

Chapter 6 provides some conclusions of this research as well as the suggestions

for future work.















CHAPTER 2
MULTIPLE VEHICLES DYNAMIC MODEL

2.1 Introduction

In order to effectively analyze orbit maintenance and re-configuration control

strategies for the satellite formation systems, mathematical models of the relative motions

between satellites are required. In the past two decades, numerous models have appeared

in the literature, with essentially the same assumptions. Typically, a leader/follower

architecture is chosen and the leader's orbit information is known [2, 21, 40, 41]. These

models can be roughly classified as shown in Figure 2.1.



Relative Dynamics Model for the Formation Flying System



PointMass Case I Rigid BodyCase



Circular Reference | Elliptic Reference | Arbitrary Reference



Clohessey Nonlinear Generalized Coupled
Wiltshire Model Attitude & Orbit


Figure 2. 1 Classification of the relative dynamics model




The most widely used model for the relative dynamics motion in the formation

flying system is the classical Clohessey-Wiltshire (C-W) equation [19-24, 26-30, 42].

Vassar and Sherwood [21], Redding and Adams [22], and Sedwich et al. [23] were the









first to apply this linear model to formation flight systems. Taking advantage of the

uncoupling between the in-plane and out-of-plane motions, a closed form controller was

designed. In 1989, Yan et al. [26] made some modifications to this linear model. One is

discretizing the continuous model, and the other is using Av (impulse) instead of normal

forces. Through these changes, the C-W solution is represented as xl = x, +FAv,

where i and i+1 represent discrete time steps, D is the state transition matrix, and F is

the input matrix. Nevertheless, the results obtained from this modified C-W model are

not acceptable for precision formation flight. For example, Vassar and Sherwood [21]

stated that the precision of the relative position was on the order of +21 Km. This

precision is too large for the types of formation flight missions being proposed.

Furthermore, this model is limited to the cases where the leader satellite is in a circular

orbit and the disturbance varies slowly with respect to time.

Queiroz et al. and Kapila et al. [31, 32, 43] used a nonlinear system as the relative

position model, which is also restricted to the condition that the leader satellite's orbit is a

circular orbit. Also imposed are the conditions that the mass of the satellite and the

disturbance experienced by the satellite vary slowly with respect to time. Lyapunov-typed

controllers are designed for these nonlinear models and the control precision is on the

order of 1 m.

Inalhan, Tillerson, and How [44, 45] in January 2002 described the relative

position motion for the case of an elliptic reference orbit, in place of the circular

reference orbit. In the noncircular case, the time dependence was replaced with the true

anomaly dependence. Inalhan's paper focuses on developing T-periodic solutions, which

refer to having the vehicles return to the initial relative states at the end of one orbital









period. Thus, the dynamic model is not general enough to study the formation system.

First it is not applicable to the singular orbit configurations. Singularity occurs when

either the inclination or the eccentricity of the leader or follower's orbit becomes zero.

Second the method only guarantees the two boundary points to satisfy the given

condition. Another approach to evaluate elliptic orbit reference case was provided by

Schiff et al. [46]. This approach, however, is limited to nearly circular orbits since it is

based on the C-W equation, which is implemented in AI Solution, Inc.'s software

FreeFlyer. This model fails for orbits with any appreciable eccentricity.

All of the studies mentioned above consider the vehicles as point masses, and as

such, only consider the formation from a relative position point of view. Ahmed et al.

[47], Wu and Chen [48], Kim et al. [49], and Xing and Parvez [50] studied the rigid body

cases and represented the attitude model using Joa + ) x Jao = M The desired state of

the attitude system is expressed as a function of the angular velocity. The basic idea is to

maintain or track a desired angular velocity profile. Based on this kind of model, different

control methods were proposed including adaptive asymptotic tracking, H,2/H, and

sliding mode control. To date, there have been no publications demonstrating a coupled

attitude/orbit problem.

Satellite Tool Kit (STK) from Analytical Graphics, Inc., provides the commercial

analytical tool for more than 30,000 aerospace, defense, and intelligence professionals. It

gives very good vision and open-loop simulation environments for spacecraft-related

research and design. However, the closed loop or control input to the spacecraft cannot be

implemented up to the most recent version. Therefore, the coupled dynamics model is

developed and addressed in this dissertation.









In Section 2.2 of this chapter, the N-body problem is simplified to the 2-body

dynamics model. Then the generalized relative dynamics model, both point mass case

and rigid body case, are derived in Section 2.3. The special cases, including all the

existing models in the literature, are also described in Section 2.3. Finally, perturbation

and their effects on the classical orbital elements are discussed in Section 2.4.



2.2 Discussion Of Multiple Vehicles Dynamic Model


Satellites experiencing formation flight constitute an N-body orbital mechanics

problem. Figure 2.2 shows a typical formation system including N satellites orbiting a

central body (Earth). The mass of satellite i is denoted by m, and that of the Earth is

M The ith satellite's center of mass is located relative to the Earth's center of mass by

vector R, and relative to jth satellite's center of mass by vector r, (r =R,-R).

Individually, the satellites are modeled as spherically homogeneous bodies with the mass

acting through their respective mass center (c.m.); however, the Earth is modeled as a

distributed mass body. Thus, for each satellite i, the gravitational potential with respect

to the Earth is given by Eq. (2.1) [51].




U= -{- [Z() kPk (cs )+k (Ckn COS m + Skm+ sin mS)Pkm (cos ()]} (2.1)
R, k=2 k=2 m=1 R,

Here R, is the magnitude of the satellite position vector R, R, is the equatorial radius of

the Earth, and (D and ( are the latitude and longitude of the satellite. Pk is the kth order

Legendre polynomial function, and Pk is the associated Legendre functions of degree n









with order m. The gravitational parameter /u is 3.986x105km3 /s2,J are the zonal

harmonic coefficients, and Ck and Skn are the sectional harmonic coefficients for

k = m and the tesseral harmonic coefficients for k m constant coefficients.



-*)1













Figure 2.2 N bodies system



The first summation in this potential function is related to the latitude variation in

the mass distribution of the Earth (see Figure 2.3a), and the second term is related to the

longitude variation (see Figure 2.4b). The darker shaded regions represent areas of

greater mass density than the lighter shaded regions. According to Prussing [52], due to

the Earth's spin, the longitudinal variation in the Earth's mass distribution can be

neglected for all satellites except geostationary satellites. Then the gravitational potential

of any arbitrary satellite i with respect to the Earth, is reduced to Eq. (2.2). Thus, for

satellites in formation systems (which are non-geostationary), the tesseral and sectoral

variations are neglected and the Earth is assumed axially symmetric; that is, the

gravitational potential is










U= -'[ I (R k) kk(COsI)] (2.2)
,R -k-.2 R.










a Tesseral b Sectoral
Figure 2.3 Zones of variations

It is well known [52, 53] that the contributions of the zonal terms (J2, J3,...) to the

gravitational potential are on the order of 103 g or smaller. Although they are small, they

have significant influence on the precision of the orbit, which is a primary interest to

formation flight. However, these contributions will be considered as perturbing

acceleration ad and will be developed in Section 2.4. What follows is the "classical"

Keperian orbit analysis.


u U= J2 R +a (2.3)


Neglecting the Earth's zonal harmonic effects and considering inter-satellite

Gm
attraction, the motion of the ith satellite is governed by Eq. (2.4), where Gm r is the


force exerted by th satellite in the formation, and ad is the disturbing acceleration

which consists of the zonal harmonic effects and other acceleration (e.g., from applied

force).









-u N Gm
R =- R, +y r +a,
S J1 9
(2.4)
_-R, N / a





Typically, formations consist of micro-satellites whose masses mn and relative

distance, are the on the order of 100 Kg and 10 Km. For a typical low Earth orbit


scenario (R &7000 Km), = Gm = 6.6722 x102km3 / s2 and =0.0081km3/s2.
r2 r2


For this typical LEO scenario, it is obvious that the contribution of the inter-satellite

attractions are even smaller than those non-spherical mass distribution of the Earth, and

so their contributions are also considered a disturbance and included in a Thus, the N-

body satellite formation problem can be simplified to an individual Earth-satellite 2-body

relative motion problem as shown in Eq. (2.5). Here 9ad includes the normalized control

efforts, disturbance and N-body attractions.


R + =R = ad (2.5)


2.3 General 2-Body Relative Motion Models


In this section, the general 2-body dynamic models for the relative motion are

derived for the formation flying system, which assumes a leader/follower or master/slave

architecture. In the development that follows, each satellite is assumed to experience 2-

body relative motion with respect to the Earth and there are no assumptions regarding to

the eccentricity of the leader's orbit.









2.3.1 Relative Motion Dynamic Model (Point Mass Case)

Consider the 2-body relative motion model shown in Figure 2.4. RL is the

position vector from the Earth's center of mass to the leader's center of mass, R, is the

position vector from the Earth's center of mass to the follower's center of mass, and p

(in Section 2.2, r, is used) is the relative position vector from the leader to the follower

satellite. An Earth-Centered-Inertial (ECI) coordinate frame is attached to the Earth

center as shown in Figure 2. 5. The angles Q, i, o), and f are the right ascension,

inclination, argument of periapsis, and true anomaly, respectively, of the satellite.

Furthermore, it is customary to define 6a = Q+a as the longitude of periapsis in the

singular cases, and u, = ac + f as the argument of latitude at epoch. Also, h and v are

the specific angular momentum and satellite's velocity, respectively. In Figure 2.6, local

vertical local horizontal (LVLH) coordinate frame is attached to the leader as shown with

x in the radial direction, z is parallel to specific angular momentum (i.e., normal to the

orbit plane), and 5 is established by the right-hand rule.


LEADER


/ \\


EARM //
( \ FOLLOWER


Figure 2. 4 Relative dynamic model






19


z






Reference Plane i -/ f Peariapis









Figure 2. 5 Coordinate system of the relative position model.





LEADER



A^l LEADER'S
)ORBE I

EA=TH

Figure 2. 6 Coordinate system of the leader satellite's frame



As defined in Figure 2.5, the position of the follower's c.m. relative to the leader's

c.m. is given by

P= RF -L (2.6)

As discussed above, each vehicle's motion is defined by a 2-body model and is

shown in Eqs. (2.7) and (2.8). Note that, here symbol q is replaced by u for the purpose









of representing all kinds of external input accelerations, includes the disturbing

acceleration due to un-modeled effects and the mass normalized control actions.


R+F -- F=u (2.7)
RF

R, +u RLL =U, (2.8)
RL

Thus, the relation motion of the follower with respect to the leader is governed by


p=R RL = [RL -( )3RFI+(F -L) (2.9)
RL RFR

Now, it is customary to coordinatize p in the local frame of the leader. Thus, the


second time derivative of p with respect to the ECI frame expressed in the LVLH is

L 00 L
L -- L L L L L L
Lp p+ x p+L+2 L x p+L iLX( LCx Lp) (2.10)


where () implies the derivative with respect to the ECI frame, and ( ) denotes the

derivative with respect to a LVLH frame. The left superscript in the vector represents the

frame in which the vector is coordinatized, and the right subscript denotes the quantity of

interest (i.e., the leader or a follower). For example, L and La) are the angular


velocity and angular acceleration of the LVLH frame with respect to the ECI frame,

L o Lo.
expressed in the LVLH frame. Furthermore, p and p are the local relative velocity and

acceleration expressed in the leader's frame. Substitution Eq.(2.10) into Eq. (2.9) yields

the governing equation for the relative position of the follower with respect to the leader.









L 00L L L X
P L 6. L X L L L _\ L) k LR n 73 i ( ) ](Ll ? \ '11\
_P=- a)x P-2 _x p- -4x( x p)+ [ R-(- ) F]+( F- L)(2.ll)
RL F

Equation (2.11) is the vector form of the generalized relative motion model for the

point mass case. Now, to evaluate this expression, the quantities defined in the expression

must be evaluated. For example, the angular accelerations of the leader and follower

L oo T
satellites, coordinatized in the LVLH frame, are 0^L = 0 0 f and



lF = 0 0 fF where fL and fF are the true anomaly accelerations of leader and


follower, respectively. The position vector of the follower and leader satellites are

defined to be FRF =[RF 0 0] and LR =[RL 0 0] To express R, in leader's

frame, the transform matrix between the two frames must be determined. Using Figure

2.5, the position vector of the follower coordinatized in LVLH can be written as Eq.

(2.12) and the rotation matrix from the follower's coordinate system to the leader's

coordinate system can be written as shown in Eq. (2.13), where C represents the


fundamental rotation matrix about the ith axis as in Eqs. (2.14) and (2.15).

F = LF RF (2.12)


C =CT C =C (lu,)C(iL)C3(QL)C 3(-F)C (-iF)C 3(-uoF) (2.13)
=LF =ECI/L=ECI/F =3 F 3

1 0 0
1C(x) = 0 cosx sinx (2.14)
0 -sinx cosx










cos x sin x 0
C (x) = -sinx cosx 0 (2.15)
0 0 1

The remaining terms of Eq. (2.11) are

L P=P, Py P (2.16)




Py fL
L LL
-L OL x p= -px L (2.17)







-L X(LLX)= f (2.18)
0
9L t)L X L P) f 2 (2.18)







-2L x p= -2p, f, (2.19)
0



If the ith element of the 1st column of CL is denoted by b then substitution of

Eqs. (2.12) (2.19) into (2.11) yields the governing equation for relative motion.


P, =p Y+L+ f2p f p[ 2+2 +(R)3RFb ]+ux, -UL
R RF
L F

P p=- L+P 2-2 + i [- ( L)3RFb2+UF -UL (2.20)
RL =F

p I -(RL )3 RFb3 + UzF UzL
R3 RF
L F











uF yF (2.21)




S= uL (2.22)


Eq. (2.20) is the generalized dynamics model for the relative motion in formation

flight. It is well know that analytic time-dependant solutions of the classical 2-body

problem cannot be obtained unless the orbit is circular. Thus, for the closed orbit case, the

quantities f, and f, in Eq. (2.20) cannot be determined. It is also well known that the

analytic solution resulting from the classical 2-body problem is dependant on an angular

measurement, which is the true anomaly. Thus, in order to be able to evaluate f, and


f,, the governing equation expressed by Eq. (2.20) must be regularized. Eqs. (2.23) and

(2.24) represents the leader's angular velocity and acceleration as functions of the

leader's true anomaly f, semimajor aL, and eccentricity eL.


/L J/aL(1-e) /a(1-e )(1+eL cos+ )2 (2.23)
rL 2 a122 )2
L a(l eL

fL -2eL (1 + eL cosfL)3 sinfL (2.24)
a 2(1 e )

Using the chain rule to change the independent variable from time to the leader's

true anomaly as the example listed in Eq. (2.25), the most generalized model for the

relative motion in formation flight is achieved as given by Eq. (2.26) or in the vector

form Eq. (2.27) [54, 55].









dg dg dfL, A
dt dfL dt
(2.25)
dg d 2
d--(g fL)=g "L+g fL
d2t dt "f+ 'f



f 2 RL RF
1 2p R
SPy -{-Py fi 2p fL LP 3 )3RFb} (2.26)
R RR
L2
p P{-1 [-( RL )3Rb3)]}
R2 L R
JL

L p i L f L -L L X f
P 2{- U L- 6 -L Lx p-2 Lx L

fL (2.27)
L L L L p RL )3 R]}
--LL x(cL xp)+ 3- -RL ) 3RF
R RF

In Eq. (2.26), the terms RFb i e[1,2,3] can be calculated from the leader's

position. This term is the ith coordinates of the follower's position expressed in LVLH.

Therefore, RFb = RL +p, RFb = py, and RFb, = pZ.


For the purpose of evaluating fL, fL, and the desired attitude as shown in

Section 2.3.5, the follower's six orbital elements: [a e Q i c f] i.e., semimajor,

eccentricity, right ascension, inclination, argument of periapsis, and true anomaly (as

shown in Figure 2.7) are need to be obtained. Therefore the relations between the position

vector R, absolute velocity vector v and the six orbital elements are needed.










Orbital Elements


:7
I,. r


-4-
~-


.Z -..'-+ u"" I
I

_j '- ,J E l.,U :
,hI- I I :,_ -. ,-


Figure 2. 7 Orbital elements



In terms of the orbital elements, the specific angular momentum magnitude is


h = jua(1-e2), the argument of latitude at epoch is u0 = f + and the magnitude of



the position vector is R =a(1- From these quantities, it is customary to write the
1+ecosf

position and absolute velocity as


cos Q cos uo sin Q sin u0 cosi
I R = R sin Q cos uo + cos Q sin uo cosi (2.28)

Lsin u sin i


-[cos f(sin uo + e sin o) + sin f(cos uo + e cos )) cos i]
vC -[sin f(sin u0 + e sin wo)- cos f(cos u0 + e cos o) cosi] (2.29)

(cos uo + e cos o) sini


In order to evaluate Eqs. (2.28) and (2.29), knowledge of the orbital elements is

required. They can be computed if the position and velocity vectors are known at some

epoch condition. The six orbit elements can be found as described below. First, the semi-


F ..,,. I .,- i .I I- I i
' .: : '- ni n 'UII LI E q u a t o r i a l
H- ,- ,. t Ecliptic

'-. ..


.,/









major axis can be determined from the position and velocity; the result is given in Eq.

(2.30). Next, from the specific angular momentum, h = R x v, the eccentricity vector and

inclination are determined from Eq. (2.31) and Eq. (2.32), where hZ is the z coordinate

of the angular momentum in the Earth frame.

2 v2
semi-major axis a= ( (2.30)
R /

1 p
eccentricity e= [vx h- R] (2.31)
,u R


inclination i = cos' h (2.32)
h

The other three elements have singularity cases.

Case 1: If i= 0, the angles Q, o), and f are not defined (see Figure 2.8). As

such, only the true longitude L = f + + is calculated here. This case includes both

e=0 and e0 .

A
1Rx
cos1 R > 0
true longitude L =R (2.33)
R-x
2r -cos R. R <0
R


Case 2: If i # 0 and e = 0, the angles ca and f are not defined (see Figure 2. 7).


[0, 0, 1] x h
The nodal vector n is calculated by n = The right ascension and argument
h


of latitude are obtained by Eqs. (2.34) and (2.35).










cos 1 77 _> 0
right ascension Q = I (2.34)
2i cos' -" < 0
n


cos1R R.v>O
R.n
true longitude at epoch u, = (2.35)
R-n
2 -cos R-v


Case 3: when i # 0 and e # 0, the right ascension Q, argument of periapsis c)

and true anomaly f are all defined and calculated from Eqs. (2.36) (2.39).


[0, 0, 1]r xh
nodal vector n [ (2.36)
h



cos-1 n 0
right ascension = n (2.37)
12, cos ny < 0
n'e

cos -- e > 0
argument of periapsis a ne (2.38)
2 cos1 q e <0
ne


cos 1- R.v>0
true anomaly f = R(2.39)
2r cos'- R v<0
R-e

To here, the regularized dynamics model of the relative position for the formation

system is derived. Considering the singularity problem, this model can be used under

both circular and elliptic reference orbit cases. In the next Section, simulation results will

show the correctness of the model.









2.3.2 Simulation Results and Discussion (Point Mass Case)

The following simulations are used to validate the generalized dynamic model

developed in this dissertation. Responses from the regularized model are generated and

compared with the Keplerian model solution as shown in Eq. (2.28).

2.3.2.1 Initial condition

Time dependant model listed in Eq. (2.20):

Initial relative position and relative velocity are calculated from Eqs. (2.40) and

F o
(2.41) below, where the subscript 0 represents the initial condition. Here RFO is the

follower's velocity with respect to follower's frame expressed in follower's frame, and

L F denotes the follower frame's angular velocity with respect to the leader's frame.

S CT C = C(UoLo)C (iLo)C(-Lo)C3 (-Fon)C1 (-o )C (-11UoFo)
=FLO =ECI/L=ECI/F =3 3

LP =F0FRF L RL (2.40)

Lo Fo L
P= RF o LLF) (CF F RF )- RL (2.41)

Regularized model listed in Eq. (2.27).

Given a set of the initial condition for the six orbital elements

[Qp,io,,co,f,,ao,eo]T Eqs.(2.28) and (2.29) are used to compute the initial absolute

position and velocity R, and v0. Through the rotation matrices, the localized velocity

and position with respect to time are calculated. The initial condition for relative velocity

with respect to the leader's true anomaly is found from the chain rule in Eq. (2.25).

2.3.2.2 Results and discussion

Due to the singularity cases discussed in the previous section, three scenarios are

simulated. The first scenario involves a zero inclination and a non-zero eccentricity of the










follower's orbit, the second scenario involves a nonzero inclination and a zero

eccentricity of the follower's orbit, and the last scenario involves both nonzero

eccentricity and inclination of the follower's orbit without loss of generalization. In all

the simulations, right ascension Q is selected to be zero, and the argument of periapsis

co is 90.

Table 2-1 shows all the cases that were tested. These cases include all

combinations of the singularity scenarios. It should be noted that only the generalized

relative position dynamics models are checked in these simulations. The dimension of the

Earth has not been considered. Therefore, the satellite may go into the Earth, which is not

practical. The derived generalized models are compared with the ideal Kepelerian model

by the norm2 difference percentage | R -R / R
-KEPLERIAN -REGULARIZED ->KEPLERIAN


Table 2-


Different cases been simulated.


Case aL (km) aF (km) iL (deg) iF (deg) fL (deg) fF (deg) eL eF
1 6878 6878 0 5 270 265 0 0
2 6878 6878 0 5 270 265 0 0.5
3 6878 6878 0 5 270 265
4 6878 6878 0 5 270 265 0.5 0.75
5 6878 6878 0 5 265 270 0 0
6 6878 6878 0 5 265 270 0 0.5
7 6878 6878 0 5 265 270
8 6878 6878 0 5 265 270 0.5 0.75
9 6878 6878 5 0 270 265 0 0
10 6878 6878 5 0 270 265 0 0.5
11 6878 6878 5 0 270 265
12 6878 6878 5 0 270 265 0.5 0.75
13 6878 6878 5 0 265 270 0 0
14 6878 6878 5 0 265 270 0 0.5
15 6878 6878 5 0 265 270
16 6878 6878 5 0 265 270 0.5 0.75
1- 6R8 68,8 0 5 270 265 0 0
IS S'S ,SON'S 0 5 270 265 0 0.5
1' ,S'S ,',SO 0 5 270 265
2W ,'Ns ,SN'xS 0 5 270 265 0.5 0.75
21 ,ts's ,xisN 0 5 265 270 0 0
22 ,NS' 5N'S 0 5 265 270 0 0.5
2-O ,,S' 0 5 265 270
24 W'ON ,,'x 0 5 265 270 0.5 0.75











CaSc I kin) a, ilmi iL (deg) iF (deg) fL (deg) fF (deg) eL eF
25 S'S s,.'s 5 0 270 265 0 0
2i I.s's .s. 5 0 270 265 0 0.5
'' S' 5 0 270 265
2s ss .s.s 5 0 270 265 0.5 0.75
2' s I5>f 5 0 265 270 0 0
I 5 i.S'S .S. 5 0 265 270 0 0.5
1 S' 5 0 265 270
S iS'S .S.. 5 0 265 270 0.5 0.75
33 6898 6878 0 5 270 265 0 0
34 6898 6878 0 5 270 265 0 0.5
35 6898 6878 0 5 270 265
36 6898 6878 0 5 270 265 0.5 0.75
37 6898 6878 0 5 265 270 0 0
38 6898 6878 0 5 265 270 0 0.5
39 6898 6878 0 5 265 270
40 6898 6878 0 5 265 270 0.5 0.75
41 6898 6878 5 0 270 265 0 0
42 6898 6878 5 0 270 265 0 0.5
43 6898 6878 5 0 270 265
44 6898 6878 5 0 270 265 0.5 0.75
45 6898 6878 5 0 265 270 0 0
46 6898 6878 5 0 265 270 0 0.5
47 6898 6878 5 0 265 270
48 6898 6878 5 0 265 270 0.5 0.75


Simulation


result 1. (Case 1 in Tabel 2-


In this case, the eccentricities of the leader and follower orbit are zero, but the

inclination of the follower is 5 (fL fF = 5"), and all the other orbit information can be


found in Table 2-1. Figure 2.8 depicits the orbits of the leader and follower satellites. In

this figure, the follower satellite is in an inclined orbit and 50 behind the leader. Figures

2.9 through Figure 2.11 show the relative position expressed in the LVLH frame. It is

seen the z direction period trend in Figure 2.8, where in half period, the follower is

above the leader, and in the other half, below the leader. For the follower is 50 behind the

leader, in LVLH coordinates, x and y directions, the follower is always later and lower

than the leader with periodic variations. Also the 3D relative position is shown in Figure

2.12. The 2-norm difference percentage shown in the results (Figure 2.13) is on the order


Table 2-1


Continued













of 0(10 16) that can be regarded as the MATLAB precision in numerical calculation.


Therefore, for the circular orbit case, the regularized model gives results that are


compatible with theppgBl& i#r model.


Follower Relative Position in x axis
Follower ^


Cl (km)


05
x axis in ECI (km)4
x 10


Figure 2. 8 Orbital information


Relative Position in y axis


0 2000 4000 6000 8000 10000 12000
time second (2 period)


Figure 2. 9 Relative position in x direction


Relative Position in z axis


0 2000 4000 6000 8000 10000 12000
time second (2 period)


Figure 2. 10 Relative position in y direction


time second (2 period)


Figure 2. 11 Relative position in z direction














Relative position


|ooo /-' i i---,---. -1
-500
~500,--
L-




-,000

-580
-20
y axis in Loca Prame ( 40 30
S-50
-620 -60 x axs in Local Frame (kim)



Figure 2. 1 Relative position in 3D


x 10-16 Difference: Kepelerian & Regularized
4


a
1 3 -


2



0
a

o
El
o
z

0


2000 4000 6000 8000
Time Seconds (sec)


Figure 2. 2 Model difference


Simulation Result 2, (Case 9 in Table 2.1):


The eccentricity and inclination of the follower satellite are both zero. As shown


in Figure 2.14, the regularized model has the difference on the order of 0(1016) The



relative position results are not demonstrated here because they are similar to the


previous case.


x 10 1 Difference Kepeleran & Regularized
4



c3

c2
(D--


(N
E


ri


0 2000 4000 6000 8000 10C
Time Seconds (sec)


Figure 2. 3 Case 9 model difference


00


10000


Simulation Result 3, (Case 20 in Table 2.1):









The eccentricity and inclination are all nonzero for the follower and leader's

orbits, and the model difference is shown in Figure 2.15. The difference percentage is on

the order of 0(10 12) This relatively larger difference in comparison to the two previous

cases is due to the additional computations required for th elliptic orbit of the follower.

Unlike the circular case, the true anomaly of the follower in the elliptic orbit cannot be

obtained from the leader's frame directly. In calculating the follower's true anomaly for

the Keplerian orbit, first the leader's true anomaly is determined from the simulation step.

Next, the regularized model uses the leader's true anomaly to calculate the follower's

absolute position and velocity. Finally, the followers' true anomaly is determined in Eqs.

(2.42) and (2.43). Computations have been conducted, which shows 0(10 12) that the

difference is coming from those above numerical steps. Furthermore, the spikes in the

plot are associated with the true anomalies of 0, .z, and 2r etc. Therefore, the

regularized dynamics model for the elliptic orbit case is also compatible with the

Kepelerian model.

RF = R + p (2.42)

(ECI)RF
cos(fF +F + + F) (c = (2.43)
RF

From the above cases, the regularized relative dynamics model is correct for the

formation flying system. The 48 cases investigated in Table 2-1 indicate that the

regularized model produce results are numerically exact to those of the time dependant

model, and thus the model appears correct. In the next condition, the special cases of the

relative dynamics model are investigated.










S102 Difference Kepelerlan & Regularized
14
S1 2


08
06

0
Z 02
0
0 2000 4000 6000 8000 10000
Time Seconds (sec)

Figure 2. 15 Case 20 model difference




2.3.3 Specialized Relative Motion Model (Point Mass Case)

Generalized dynamic models of the relative position motion in the formation

system have been developed in Section 2.3.1. In this Section, three specialized models

are studied. The first two have been heavily used in applications and literatures which are

the C-W equation and the in-plane nonlinear circular reference orbit equation. The third

one is the in-plane nonlinear elliptic reference orbit equation. Here in-plane represents the

case where the leader and follower satellites are in the same plane.

2.3.3.1 Linear model (Clohessy-Wiltshire equation)

If the reference orbit is a circular orbit, the general 2-body relative motion

dynamic model derived in Section 2.3.1 can be simplified to a linear system, which is the

well-known Clohessy-Wiltshire (C-W) equation. According to Section 2.3.1, the second


derivative of the relative position with respect to the ECI frame can be expressed in a

LVLH frame, here, i.e., the leader's frame as Eq. (2.44). For the general 2-body motion

/1R
R+ RF R =F, let g(R) = RR, then Eq. (2.45) can be achieved. Substitution of the



relative position p = R1 -RL into Eq. (2.45) and using a Taylor series expansion about









reference position R,, Eq. (2.46) can be derived. Neglecting the higher order terms and

utilizing R, = g(R,)+ u, the linear model is derived as Eq. (2.47). Now the equation is

written in the reference frame, i.e., leader's frame, using Eq. (2.48) and UFL = UF -L

L .. L L
L -- L L L L L
p= p,+ p+Lx Lp +2L cLx pL+ Lx( LLx Lp) (2.44)


RF = g(RF) +u (2.45)

g(R,) 2
p + LR = g(RL)+ p + O(p )+ F (2.46)


= (RL P+uF -L (2.47)
RL

L L L L L L L g(RL)

RL T


In the circular reference frame case, oL = 0 fL = [0 0 0],


200
L L 0 L =[RL 0 1 and Og ) /
L 0 0 j =[, O O] and 0 -1 0 Finally, the C-
-L0 0 -1

W equation is derived as Eq. (2.49), where matrices A and B are shown in Eqs. (2.50)

and (2.51).










P PF Px uFLx
py =A pY +B py + uFLY (2.49)
= _P_ P. MFLz


3fL 0 0
A 0 0 0 (2.50)
S2
0 0 -fL


0 2fL 0

B= -2fL 0 0 (2.51)
0 0 0


The C-W equation was derived in the nineteenth century [52], which originally

described the motion of the Moon relative to the Earth. Although it is different from the

formation flying application, they both consider small displacements relative to a known

reference motion. This linear model has been used extensively in formation flight models

[19, 20, 21, 22, 23, 24, 26, 27, 28, 29, 30].

2.3.3.2 Circular reference orbit nonlinear model

Due to the uncoupling between the in-plane and the out-of-plane motion,

nonlinear co-planar relative position equations with respect to a circular reference orbit

will be shown in this part. In this case, the leader satellite is limited in the circular orbit.

Therefore, proper coordinates are set to let the inclinations for the leader and the follower

satellites to be the same, i.e., i = iL = 0. Here the right ascension is assumed to be

zero. Based on these conditions, from Eq. (2.20), the circular reference orbit in-plane

model can be simplified to Eq. (2.52), which has the same results as derived from the

Keplerian model.












P = Px fl+2p L [RL- (R )3 RF (coso cos MOF + Sin L Sin UF )] + UF xL
L F (2.52)

Py = Py f 2p3 R1 (( RL )3R (- sin uOL cos uOF + cos sin )]+u yF yL
I R L RF

2.3.3.3 Elliptic reference orbit nonlinear model

Just as shown in section 2.3.1, in elliptic reference orbit case, it is not convenient

to represent the model as a function of time, because the angular velocity of the leader

satellite is not a constant. Therefore, the true anomaly is used as the independent variable.

Based on the dynamic model derived in the previous section and using the chain rule

dg dg df the governing equations for the in-plane elliptic reference orbit relative
dt df dt'

motion can be re-written as Eq. (2.53), where the notations are listed in Eqs. (2.54)-

(2.57). As Section 2.3.3.2, this model is also an in-plane special case derived from the

general relative dynamics model.


px = 7 {-fLpX + 2(f,)2 p + [(fL)2 /- ]px +fLp+ cRLx + UxF L}
(f-) (d +d2)32
L)d + 2)/ (2.53)
1 +
Py = {-- ^ 2( f)2 P + [( 3/2 LPy +C2Ry +UF UyL}
(fL)2 l (d, +d)

2/ 7?[2 [aL e(1 Le )]2
(1 + e, cosf)2

aL(1 -e)
d2 =p p2 +p 2 (2.55)
(1+ e cos f)









d3 =hL =u(1-e ) (2.56)

c, =c= -2 [(d, +d2)3/2 dJ2] (2.57)
S2 Jd (d + d2)3/2

Up to now, three special cases of the relative motion models have been

introduced. In the next section, how precisely the C-W equation is are shown by

simulation.

2.3.4 Simulation Results and Discussion (Point Mass Case)

Simulation results show the difference between C-W equations and the

generalized relative position model. Results demonstrate the necessary of using nonlinear

model instead of the C-W linear model in the formation flight system. The in-plane

circular and elliptic reference orbit model are not compared, since it is identical to the

generalized relative position model under the in-plane assumption. The leader is in an

equatorial orbit with 6878 km radius and 0 true longitude. The follower is also in the

equatorial orbit, with 0 true longitude. But the follower's radius is 6888 km. Initial

relative velocity with respect to the leader's frame is calculated by


Vf V R= R


Figure 2.16 and Figure 2.17 show the relative position results from the C-W

equation and the generalized relative dynamics model. In Figure 2.18, the difference

between two models is shown. The difference is calculated by the 2-norm between these

two models at each time second. C-W equation has analytic solution as shown in Eq.

(2.58), where n = fL. p, and p,0 are the initial relative position and velocity for

i e [x, y, z]. Although C-W equation is easier to calculate, it can only be used in the













circular reference orbit. Furthermore, due to the truncation of the high order terms, the


solution is not periodic, which can be seen from Eq. (2.58). These are the reasons why the


nonlinear model is adopted in this dissertation.



Px = PA cos nt + P" sin nt
n


p = (4 3cosnt)po + AO sinnt + 2(1- cosnt)P,, (2.58)
n n

2 4 sin nt 3nt
p, = 6(sin nt nt)po + Py (1- cos nt)p.o + PYo
n n


Comparison of CW and Generalized Model Comparison of CW and Generalized Model
60 100
CW
0 Generalized Model
S100
I .2 -200
o 20200
.: -300
0 -400
W -500
S-20 Generalized Model -0
600
-700
-40 -
-800
-60 -900
0 05 1 15 2 0 05 1 15 2
Time Seconds (3 Periods) x 104 Time Seconds (3 Periods) x 104


Figure 2. 16 Relative position in x direction Figure 2. 17 Relative position in y
direction




Difference Between CW and Generalized Model
120

100

80

60

40

20

0
0 0 5 1 1 5 2
Time Seconds (3 Periods) x 10


Figure 2. 18 Difference between CW and generalized model









2.3.5 Relative Motion Dynamic Model (Rigid Body Case)

In section 2.3.1, the satellites are regarded as point mass vehicles and the relative

position model is derived. In this section, the rigid body case is studied, which includes

both position (orbit) and orientation (attitude) motions.

To date, most literatures [19, 20, 21, 22, 47, 48, 49] regard the orbit and

orientation of the formation flight as an uncoupled system. The relative position and

attitude motions are studied independently. However, this is not necessary true. For

example, missions as EO-1 and LandSat7, which were required to collectively point to a

single spot on the surface of the Earth, once or twice a day. Furthermore, one of the final

objectives of the formation flight system is pointing and positioning managed

collectively. Therefore, the coupled position and attitude models must be studied.

In this dissertation, the objective of pointing managed collectively is simplified to

be the following. All the satellites in the formation flight system are required to point to

the nadir of the leader satellite. Also, all the satellites are assumed to be gravity gradient

stabilized and thus are pointing to their own nadir points under the passively stabilized

condition. Based on the above simplification and assumption, only the desired

orientations of the follower satellites are coupled with the real time relative position.

In this section, first, the desired attitude model is derived. Then, the attitude

governing equation is studied. For the convenience of using sliding mode controllers,

which will be discussed in Chapter 3, the Euler angle representation for attitude are

adopted.

2.3.5.1 Desired attitude

For a gravity gradient stabilized satellite, the body frame is chosen as shown in

Figure 2.19. For instance, shown in the follower's body frame, axis z points to the









follower's nadir point NF, k is perpendicular to axis and in the follower's orbit plane,

and ) is determined by the right hand rule. The desired attitude of this follower satellite

is that the z direction points to the nadir of the leader satellite NL. Define two error

angles, an elevation error EL along latitude and an azimuth error AZ along longitude. In

order to calculate these error angles from the real-time relative position, three theories

coming from the spherical geometry are used.


z




















Theorem 1 [Spherical Pythagorean Theorem Ref. 56] In a right triangle AABC

on a sphere of radius R with a right angle at vertex C as shown in Figure 2.20, angles a ,

/Y and y are facing sides BC, AC and AB separately. The lengths of sides BC, AC
and AB are a, b and c, thenOLL
cos( ) c ) c ( )LADER























R R R
Figure 2. 19 Satellite attitude body frame



Theorem 1 [Spherical Pythagorean Theorem Ref. 56] In a right triangle AABC

on a sphere of radius R with a right angle at vertex C as shown in Figure 2.20, angles a ,

p and /are facing sides BC, AC and AB separately. The lengths of sides BC, AC

and AB are a, b and c, then


cos( = cosa -cos(-) (2.59)
R R R









Theorem 2 [Spherical Sine Theorem, Ref. 56] Let AABC be a spherical

triangle on a sphere of radius R. Let a, b, and c denote the lengths of the sides, and

let ZA, ZB and ZC denote the interior angles at each vertex. Then


sin( ) sin( ) sin( )
R R R (2.60)
sin(ZA) sin(ZB) sin(ZC)

B




zqz/


A b
Figure 2. 20 Spherical right triangle




Theorem 3 [Ref. 56] Suppose that AABC is a right triangle on the sphere of

radius R with a right angle at C. Then cosA = sin Bcos(a/R).

Figures 2.21 and 2.22 show the definitions of the error angles. Let i =iF -iL be

the inclination difference between the follower and leader orbits. The curve line NFA is

on the longitude line of the follower satellite and elevation error angle LF is facing this

arc line. Point A is the intersection of NFA and the leader's orbit plane. Along the

latitude line of the leader's orbit, NL is the nadir point of the leader and azimuth error

angle faces the arc line ANL In Figure 2.21, the angle between ONF and OA is L. In

Figure 2.22, the angles between OB and OA, OB and ONL are defined to be aF and

a, separately. In order to obtain the desired attitude, the follower satellite first rotates











along its longitude from the point NF to point A about its body axis i through LF.


Then, i points to NL through angle ac -aF after rotating about its body axis j In the


Figure 2. 23, the angles' definitions are summarized below. These angles are to be used


in calculating the error angles.


FOLLOWER


Y


X ORE IT
** X


Figure 2. 21 Desired attitude I


XMrI


YVBr


LEADER


Y










Figure 2. 22 Desired attitude II


Angles Definition Restriction
LF1 Angle between lines OD and ONF -iF LF < iF
LF2 Angle between lines OD and OA -iL < LF2 LF Angle between lines ONF and OA LF = LF1 -LF2
UOF Angle between lines OB and ONF 0 < uOF < 2;

aF Angle between lines OA and OB 0 < a < 2f;
aF1 Angle between lines OD and OB 0 < aF < 2l ;


FOLLOWER
ORB IT


LEADER
ORB IT


LEADER
ORBIT


Figure 2. 23 Desired attitude III


Now, the error angles are to be calculated from the real time relative position

information. When the relative position of a follower satellite (i.e., p", py and p,) is


found, the orbit elements of this follower, such as u0F, iF and QF, can be determined as

discussed in Section 2.3.1. In Figure 2.20, using the spherical sine theorem, in triangle









ABDNF, LF1 can be found through Eq. (2.61). Also, in this triangle, aF can be

calculated by spherical Pythagorean theorem through Eq. (2.62). Note that, the quadrant

of the cosine function is determined by the true anomaly uF.

LF1 = sin (sinuOF siniF) (2.61)

aF = cos (cos uOF / CosLFl) if 00 (2.62)
aF1 = 2 Cos- (cos uOF / COS LFi) if T < UOF < 2(.

Here, OF needs to be added into aF, because the follower orbit's right ascension

of the ascending node may change due to the control effort and the orbits of the leader

and follower may not cross at the same point. Therefore, let a.F2 = aF1 + QF. Based on

the theorem 3 and spherical sine theorem, in spherical triangle ADBA, LF2 and aF is

given by Eqs (2.63) and (2.64). Therefore, LF is obtained by LF = LF1 -LF2. From these

two angles LF and aF, the error angles are to be derived.

LF2 = tan(sinaF2 taniL) (2.63)

SaF- =cos (cosLFz2 COSCF2) if 00 (2.64)
acF, = 2 co o (coLF2 -COSaF2) if 7
Two rotations should be taken to achieve the desired attitude. According to the

Figure 2.24, first, follower satellite rotates EL about the body axis x, and second it

rotates AZ about the body axis 5 The following steps are used to generate the real time

EL and AZ. All vectors are expressed in ECI frame for the convenience of calculation.












FFOLLOWER
Election FK




LEADER







SX
IL









Figure 2. 24 Error angles in the coordinate system. EL is the elevator error angle,
whereas, AZ is the azimuth error angle.



R,
(ECI) OA = R Q(-L)R (-iL)R3(-aF) 0 (2.65)

0


RL + px
(ET)FO = -R3(- ,L)R (- iL)R3 ( ) (2.66)
A

(ECJ)FA = ((ECI)FO +(ECI)OA) (2.67)

(ECI)FA (ECI)FO
EL = cos-' [ (EC)FA (EC)FO (2.68)
(ECI) FA | (ECI) FO


Considering about the rotation direction when the follower satellite is in different

quadrant, the satellite rotates about i axis in negative direction when the true longitude

of the follower is in [0,;r) and positive when [i, 2r) Now, the error angles between

actual and desired attitude have been derived as shown in Eqs. (2.69) and (2.72).









(ECI)FA (ECI)FO

(EC)FA (EC)FO
-"(EdJ) (Ec-)F,] (2.69)

EL =-cos (ECI)FA (E ] f 0<



EC)ONL = R (-L)R (-iL)R (-aL) 0(2.70)
0

(EC) FNL = (EC) FO + (EC) ONL (2.71)

(ECI) FA (ECI) FN
AZ = cos 1[ ) (2.72)
(ECIJ)FA| (ECI)FN L

After the error angles between the desired attitude and actual open loop attitude

are calculated, the desired attitude can be obtained through the direction cosine matrix

form. For the gravitational gradient stabilized satellite, the unforced attitude only has one

principal axis rotation, which is the rotation against the body axis y with angular

velocity, which is related to the orbit period. Let's denote this open loop Euler angle to be

0 Then the desired attitude direction cosine matrix is calculated in Eq. (2.73), where

subscript i e [1, 2,3] means the rotation with respect to ith body axis.

CD = 2(AZ)C (EL)C2 (O) (2.73)

2.3.5.2 Attitude model

The attitude motion of a rigid body is governed by

JOBODYDY OY BODY BODY (2.74)


where J is the centroidal inertial dyadic of the body, OBODY is the angular

velocity of the body coordinatized in the body and -BODY is the extErnally applied torque









(coordinatized in the body). In terms of Euler angles, the angular velocity can be written

as Eq. (2.75), where the rotation sequence is chosen to be 1-2-3 through

= [EL, AZ, 03] = [01,023 Rotation 1 is for the elevation error angle and rotation 2

is for the azimuth error angle. Rotation 3 is used here only a supplement in order to get

three differential equations and 03 3= 03 = The angular acceleration is derived as

shown in Eq. (2.76).

0BODY =8x1 + 822 + 83 3
cos 0 sin 3 0 0 cos 03 sin 3 0 cos 2 0 sin 2
=033+ -sin3 cosO3 0 02 + -sin3 COS03 0 0 1 0 0
0 0 10 0 0 1 -sin2 0 cos 02


= 2cos 3 + -sin 3 cos 3 0 0 (2.75)
03 0 1 -0, sin 02
02 sin O3 +01 cos2 cos 3 cosO2 cosO3 sin 03 6O 0O
= 02 cos03 0cos02 sin = s in 03 cos 80 0 o2 =S() 2
03 -8 sin0 2 -sinO2 0 1 03 L03]



cos02 cos 3 sin 03 0 01 0203 cos 3 -102 sin 02 cos 3 -0103 cos02 sin 03
W BOD = cos 2 sin 03 cos 83 0 2 + -02 sin 03 + 82 sin 02 sinO3 0- O cos02 cosO3
-sin02 0 1 3 -0102 cos 02
=s() + S,(o, o)
(2.76)

where = is a trigonometric function of the last two Euler relations and d is a

column matrix of the Euler rates. Substituting Eq. (2.75) and (2.76) into Eq. (2.74) results

in










S= (JSy) [-JS, SOJSB + TBD ] (2.77)

The attitude model in Eq. (2.77) is used for the large angle case. In the small

angle perturbation, the infinitesimal angles can be used to linearize the attitude governing

equations. For any small disturbance 30 = [0, 8302 803 ], the actual angular velocity

and acceleration of the follower body can be written as O BODY = QBODYO + M and

)BODY = BODYO +08, where subscript 0 represents the stabilized point. Neglecting the

high order term, the small angle perturbation equations can be written as the following

linear system 0 = C30 + DT. Here, C and D are constant matrices determined by the

stabilization point and T is the control torque. Furthermore, because the relative position

model is in terms of true anomaly instead of time, the attitude governing equation also

needs to use the chain rule and be changes from the time dependant version to the true

anomaly dependant version.

2.3.6 Simulation Results and Discussion (Rigid Body Case)

Simulation Result 1:

The leader satellite is in an equatorial circular orbit, with a 500 km altitude and

the true longitude to be 0. The follower satellite is also in a 500 km circular orbit with a

50 inclination and a -50 true longitude initially. Assume the satellites are gravity

gradient stabilized and thus naturally points to their respective nadir. The open loop

simulation results are show in Figure 2.25. Figure 2.25 shows that the elevator error angle

is periodic. This is expected as shown by the angle definition in Figure 2.24 and the

chosen orbits. The follower satellite rotates about x in the negative direction in the first

half of the period, and then positive in the second half of the period. The azimuth error










angle is also periodic. The maximum value in the elevator error angle can be shown to be

on the order of 0.87 radians and this is confirmed in Figure 2.25.

Simulation Result 2:

In this scenario, the case when the inclination of the leader is non-zero is

evaluated. The only differences in initial condition from the simulation 1 are the

inclination of the leader is 50 and that of the follower is 100. The simulation results are

shown in Figure 2.26 and are identical to those obtained in simulation case 1. Because the

relative inclination between the leader and the follower has not been changed, thus the

results are not changed.

1 Simulation Results 1 Simulation Results 2



o a0

Thetal Thetal
SI --Theta2 Theta2
\-0 5 o 5

-1 .1
0 2000 4000 6000 8000 10000 12000 0 2000 4000 6000 8000 10000 12000
Time Seconds (2 Periods) Time Seconds (2 Periods)

Figure 2. 25 Theta 1 (EL) and Theta2 (AZ) Figure 2. 26 Theta 1 (EL) and Theta2 (AZ)
in radians in radians



Simulation Result 3:

The case when the follower is in an elliptic orbit is shown in Figure 2.28. In the

point mass case, the Earth dimension ha not been considered. Here in the rigid body case,

the Earth dimension is included. So, the elliptic orbit needs to be chosen carefully to

ensure that the satellite will not fly into the Earth. The semi-major axis of the leader and

the follower is chosen to be 13868 Km. Inclination of the leader is 0, and that of the










follower is 5. Also the eccentricity of the follower satellite is 0.5. True longitudes of the

leader and follower are 0 and 50 respectively. As shown in Figure 2.28, 0, goes to a

large angle when the leader's true longitude arrives kzr, where k =0,1,2.... In the

elliptic case, the error angles are also periodic.

Simulation Results 3


SThet
+- Thet







0 1 2 3 4
me Seconds (2 Periods o

Figure 2. 27 Thetal (EL) and Theta2 (AZ)



2.4 Perturbation

Until now, all the equations and simulations have not included disturbance

effects. Actually, a satellite in an Earth orbit experiences perturbations, although the

absolute value is small. As stated in Prussing [52], "The perturbations come from the lack

of spherical symmetry of the Earth, the attraction of the Moon, and Sun, and atmosphere

drag if the satellite is in the low orbit." The perturbation forces are normalized and

expressed as accelerations (force per unit massKN/kg) and coordinatized in LVLH

frame. In this dissertation, only perturbations associated with the atmospheric drag and

the effects of the Earth oblateness are considered.









2.4.1 Atmospheric Drag

The drag force acts in a direction opposite to the velocity of the satellite. For a

satellite is in a low Earth orbit (LEO), it experiences a drag due the Earth atmosphere

which can be represented as Eq. (2.78) [52].

1 CA 2
ard =- -p vv (2.78)
S 2 m

Here, p is the atmosphere density, CD is the drag coefficient, A is the cross-

section area of the satellite, v and v are the velocity magnitude and unit vector of the

satellite, and m is the mass of the satellite. The quantity CDA/m is called the ballistic

coefficient of the satellite. Subscript ad represents atmospheric drag. Therefore, the

1 CDAL 2
atmospheric drags for the leader and the follower are a dL =--p L L and
2 mL


a p CD F 2 F respectively. The relative atmospheric drag for the generalized
a 2 m F

motion in a formation system is calculated by Eq. (2.79), where t i e 1, 2,3 is determined

by CLF, aaF, and a dL Here, L = [eL, IeL, e LT are the leader's LVLH unit vectors. The

table below listed the coefficients for the atmospheric drag model. The altitude of the

satellite is on the order of 500 Km in the simulation. The dimension and weight of the

satellite are chosen from ST-5 [51].

Atmospheric density p Drag Coefficient CD Area A Mass m
4.89 x10 l kg/m3 (500km altitude) 2 (typical) 0.286x0.542m2 19.5kg


adF adL = tleL +t2eoL +t3e


(2.79)













Z A=


A
A







EA=TH









Figure 2. 28 Atmosphere drag



2.4.2 Earth Oblateness

As was discussed in Section 2.2, the averaging effects of the Earth's spin results

in a gravitational potential of the satellite with respect to the Earth as expressed in Eq.

(2.80) [52]. R is the distance from the Earth center to the satellite, R, is the Earth


equatorial radius, Pk is the Legendre polynomial function of the kth order and D is the


colatitude. The first three coefficients are J2 =0.00108263, J3 =-0.00000254,


and J, = -0.00000161.


u= [1- Z ekJP(os)] (2.80)
R k=2 R


Since J2 is significantly bigger than the others, the satellite's potential with

respect to the Earth is approximately as










U= [1 (R)2 Jp ( ()] (2.81)
R R

where the Legendra polynomial is


P2 (cos ) = (3 cos2 1)
2

z
Can be determined as Figure 2.29, the quantity cos O is simply where Z and
R

R are Z component and the magnitude of the satellite's position respectively. Thus, the

potential including oblateness effects is

pU 1 Re Z
U= 1- ) J2[3 1)]} (2.82)
R 2R R2

In the expression for the potential, the first term is related to the central force,

which has been previously considered. [52] Thus, second term is the disturbance related

to the Earth's oblateness

j 2R 1 Z2
U J 2[ (3 1)] (2.83)
R 2 R


EARTH
1 -



2. 29 Te p n ss fm te Eh

Figure 2. 29 The perturbation stems from the Earth oblateness


A









The perturbating acceleration can be derived as Eq. (2.84) [52], where the

subscript obl represents the Earth's oblateness. In order to accommodate the

perturbations into the general dynamics model for the relative position in the formation

system, the difference between the follower and the leader's oblateness pertubation need

to be derived as shown in Eqs. (2.85) and (2.86). Therefore, the J2 term perturbation for

the relative position model is written as Eq. (2.87) where p, i e1,2,3 is determined by

Eqs. (2.85) and (2.86).

OUJ2 OU
Rol r + Z
r Z Z (2.84)
3-(3 15Z2)r J2 3Z (
= -pJ 2R ( 4 2R6 r 2R 5 Z


3 15Z2 3Z
qobF =-J2R2(-4 5_ rF 2e)e Z (2.85)
2RF 2RF RF

3 15Z2 3Z
qobL 2= -jR2 ( )rL pJ2R Z (2.86)
2R, 2R^ R,

aoblF QoblL = PlrL + POeL + p3zL (2.87)

Include the J2 and the atmospheric drag terms derived in Section (2.4.1), the

generalized relative position model derived in Section 2.3 becomes Eq. (2.88).










1 2 o 2 P RL
P {-p L+L2p fL+pf L xpfL+ -[RL -( )3 RFbl]+t +P
R R

1 2 RL
S = {-P, fL -2p fL -A + L +[- L + ( ) ]+2 +pt } (2.88)


fL



2.4.3 Perturbation Effects on the Generalized Dynamics Model

Simulation 1: In this section, the effects of the perturbations are shown. The

satellite's orbit is a circular orbit with 0-degree inclination and 500 km altitude. As

shown in Figure 2.31, the altitude of the satellite is decreasing due to the atmospheric

drag. Seen in Eq. (2.78), the atmospheric drag always acts on the velocity direction and

decrease the velocity due to the negative sign. From the inverse law model (or Keplerian

model), less velocity causes lower orbit and degrades the altitude of the orbit. However,

the J2 effect causes periodical variations on the satellite altitude. J2 comes from the Earth

oblateness, which is caused by the non-uniform distribution of the Earth density. But the

total mass does not change. Therefore, satellite's altitude changes periodically, but the

maximum altitude maintained (as shown in Figure 2.32). For the simulation condition,


where i=0, the J2 acceleration aobl can be simplified as --x J24 which is in the
2R4

negative radial direction. This acceleration causes the satellite to be driven inside.

Furthermore, J2 changes right ascension Q and argument of perigee as shown in Eqs.

(2.89) and (2.90). Therefore, the right ascension node and argument of perigee change

due to the J2 effects. In order to easily show the effects using simulation, J2 coefficient is












made to be 100 times larger than the correct number. Figure 2.33 shows the combined


effects from the above two reasons.


S2 = -2.06474x1014a 7/2(cosi)(1 e2)-2 (2.89)



1) 2 = 1.03237 x 1014 7/2(4- 5sin2 i)( e2)-2 (2.90)


Altitude degradation due to atmospheric drag Altitude degradation due to J2
6880
Inversesquar- Inverse square model
Inersesuareodl J2 perturbation
Atmosphenc drag p i i
6875 I
f~r ""-- 6 5 ; i i i t t


s6877 9E 6870

O I = i
,877.899 6865 .
ii
3"n9i 6860
i i! 'i !! ;i i il i ;1


,i ---- ^---,-------- 6856 ---------------
S0 1 2 3 4 5 6 7
ie beni CA:.Ji 0 5 I:pert I Time seconds (10.5 periods) x 10

Figure 2. 30 Atmospheric drag effect Figure 2. 31 J2 perturbation effect

Orbit cegradation due to J2 Perturbation
-._ Desired Orbit
-" .-- J2 Perturbed Model






0 3 --- -- "" ^ ^ ..E rrj| ^ ------''^ '. '- "
08
0 5 .-2 a- ------






2000 '- _"^ :- "-. .. ..-. --'_. "'- 04- ---0







7900 -02---- .-------_ """ 0
0 00 .--- "
06 Absolute poson y drect (km) -00 00 Absolute potl In direction (km)
n4 --s----








Absoutepostio in diecton krn) -E0C _000 _E00 -600 400 A~slut poitin inx drecion(kr


Figure 2. 32 Combined perturbation effect










Simulation 2: In this part the effects of disturbance on the relative dynamics (Eq.

(2.26)) model are shown in Figures 2.33 and 2.34. The leader satellite is in a 6878 km

equatorial orbit. The follower satellite is the same radius circular orbit as the leader,


except with 50 inclination. The true longitude of the leader and follower are 0 and 5,

respectively. The 2-norm difference between the case without perturbations and the case

with perturbations are shown. The absolute errors are very small which are on the order


of O(10-4) in the relative position model and 0(10 6) in the relative attitude model.


Since the errors are increasing with the time, these perturbations must be considered.

x 1-4 Pertutration's Effects onthe Relatie Position x 10 Petrtaatiris Effects on the Dsired QC terion
5 1.8

1.6
4 1.4

1.2
^3
E 3 1
J 0.8
-2
0.6
1 0.4

0.2

0 2000 4000 6000 8000 10000 12000 0 20 4 6000 8000 1000 120i
Tine second Tirresecord

Figure 2. 33 Error in position Figure 2. 34 Error in attitude




2.4.4 Perturbation Effects on Keplerian Elements

In the previous section, simulation results have shown the perturbation effects on

the absolute and relative dynamic models. In order to access the perturbation's effects on

the six Keplerian orbital elements, the Gaussian Form of Lagrange's Planetary Equations










is derived according to reference [57]. Here, the brief derivation and explanation are

illustrated.

The orbital elements vector is listed in Eq. (2.91) and the detailed information is


given in Table 2-2.


a = [ / i o a e M]A


(2.91)


Table 2- 2 Classical elements for an Earth satellite
Quantity Specified Symbol Definition Special Cases
by
Reference Time Epoch To Time for which elements are
specified
Size Semimajor a Half the long axis of the ellipse
axis
Shape Eccentricity e Distance for center of ellipse to 0: circular orbit
focus divided by semi-major axis;
dimensionless ratio.
Orbit plane Inclination i Angle between orbit plane and 0 :equatorial orbit;
earth's equatorial plane 90 polar orbit
90" : polar orbit
Right 0 Angle measured at the center of 0 : the ascending
ascension the earth in the equatorial plane node is at the
of the from the vernal equinox eastward vernal equinox.
ascending to the ascending node.
node
Orientation of Argument cD Angle measured at the center of
orbit in the of perigree the earth in the orbit plane form
plane the ascending node to perigee.
Location of the Mean M 3600 x At/ P where P is the 0, when the epoch
satellite in its anomaly at time is perigee.
orbit epoach period and At is the time
difference between the epoch and
the last perigee passage before
_epoch
In order to apply the Lagrange's Planetary Equations (LPE), the mean anomaly


M must be expressed as


M = nr, where n = P is the mean motion, and r is the time
Ja


of pericenter passage. The rotation matrix from the Earth frame to the local (i.e., the


satellite) frame is given by










C 1 m1 m2 m3

n1 n2 n3
cos 0 cos c sin Q sin a cos i cos 0 sin a sin Q cos a cosi sin Q sin i
= sin Q cos o + cos sina cosi -sin Qcos + cosQcosa cosi -cosQsini
sin o sini cos osin i cos i

(2.92)

Then the satellite's absolute position and velocity vectors can be written as Eq.

(2.93) and Eq. (2.94).

acosE-e
R=C bsinE (2.93)
0


a f sinE/(1 e cosE)

= C bfcosE(1-ecosE) (2.94)
0



E l-e f
where E is calculated by tan- e- tan
2 1+e 2

Assume the disturbance scale function is W. The system equations of the 2-body

dynamic motions under perturbation can be written as Eqs. (2.95) and (2.96). Let the

homogenous solution for the system equations to be Eqs. (2.97) and (2.98). Then the

basic method of variation of parameters is used here, the solution of R and v are

assumed to have the forms shown in Eqs. (2.99) and (2.100).









dR
=1= (2.95)
dt

dv + R W[]T (2.96)
dt R3 OR

R = R(t, a) (2.97)

v= v(t, a) (2.98)

R = R(t, (t)) (2.99)

v= v(t,a(t)) (2.100)

Taking the derivative of the vector R and v and substituting these derivative into

Eqs. (2.95) and (2.96), the partial differential equations are derived as shown in Eqs.

(2.103) and (2.104).


IR cR OR da
dt 8t ca dt

dv yv Ov da
dt at Oa dt

OR da
=O0
da dt

yv da OW r
Oa dt OR


(2.101)


(2.102)


(2.103)


(2.104)


Using [ ]o left multiplied with (2.103) subtracting [ ]T left multiplied with
(2.104), the LPE are derived as Eq. (2.105)
(2.104), the LPE are derived as Eq. (2.105)









Lda
Sdt


dW
[ ]'
Sa


(2.105)


where L = [ ] [ ] In the mean time, reference [57] gives the derivative of
Oa Oa Oa Oa

the position and velocity vectors with respect to the six orbital elements as Eq. (2.106)

and (2.107).


=a(1-e) 1,
80
0

OR 1
8= a(1 e) m2
1n2

OR
S-a na m
8e

nL

2
S=fb/(l-e) 12


R1
= a(1 e) sin o m,
8i
n,/

LR i 3bnt 12
S=(1-e) mn, -- m
Oa 2q



In2
= b /(1- e) m,
LiJ


1 2 -1 12
/F 8l 3f t /b
= j bcos0/(1-e) m3 3 bf m 2 (2.107)
Oi 3 a 2(1-e)2 1 2a(- e) 2
L/3 1, /2
L['2 1'

8e b(1-e) 2 A (1-e)2


Combined the LPE with above partial derivative equations, the governing

equations are derived as Eq. (2.108). Notice that the matrix L is nonsingular if e # 0 or

e and i 0. Also, according to reference [57], the disturbance scalar function's

partial derivations with respect to the six orbital elements are listed in Eq. (2.109).


(2.106)


O1V
av
= f b/(1- e) m,
-nj
da0
L.7iJ









1 OW

fab sini ai
cosi OW b OW
+-
f absini i f ae Oe
b OW b2 OW

fa3e fo) a4e 9


0
9W T 9W
- u d Rcosi =
-RcosOsini


0
OW T 0
Ud 0
di
r sin 0



F0o
OW T
aa l1d L
00)
0


W
ae





OW
- 2
aA


1 OW cos

f ab sini f abs
2 OW
na OA

2 OW b2 OW

fa 9a fa4e 9e


R 3aft .
e sin f
a 2b

S3a ft ( f)
*Ud (1+ecos f)
2b
0



-a cos f
T aR
ud (1+ )asinf

0
2
ae .
--sin f
b
2
Ta
d -(1+ecosf)
b
0


i SW

ini )


(2.108)


SJ (2.109)

Finally, Gauss form of LPE is derived to be Eq. (2.110), where p = a(- e2),


R = ,ec and h = p From these system equations, it is easily to find that the
1+ecos f

perturbation force ud =[udr udo ud ] affects the six orbital elements. Therefore, the

orbital elements are not constant under the perturbations.












0 0 R sin() + f)
hsini
0 0 R cos(o + f)
h
0
1 pcosf (p+R)sinf Rcos isin(a+f) -Ur
dr 0
d a he he hsinidr
- dt a 2ae sin f 2ap 0 UdZ 0
e h hR 0
M psinf (p+R)cos f +re 0
h h
(pcosf -2Re)b b(p + R)sin f
ahe ahe

2.5 Conclusion

The generalized dynamics model has been developed with the inclusion of

perturbation effects. The developed dynamic model has been shown to reproduce the

"special case" models currently in the literature. The derived generalized dynamics model

can be used in any closed reference orbit (i.e., circular or elliptic reference cases). It is

shown that for the scenarios considered in this dissertation, the desired attitude is coupled

with the real time relative position. Perturbation effects from the atmospheric and Earth's

oblateness are studied and their effects on the Keplerian model, relative position, and

orbital elements are evaluated. Although those effects are small in the magnitude sense,

their cumulative influence is significant and needs to be considered.














CHAPTER 3
CONTROL STRATEGIES: ORBIT

3.1 Introduction

In Chapter 2, the generalized dynamic models, including both the point mass case

and the rigid body case, have been intensively investigated. In Chapters 3, 4, and 5, the

controllers are designed in order to maintain and/or re-configure the formation. The

controllers for the formation system can be treated as a two-level architecture, in which

the higher level dealing with the systematic control and the lower level related to the

individual leader-follower pair. The higher level of the controllers will be studied in

Chapter 5 where the information about how to manage the whole formation system is

discussed. In Chapters 3 and 4, the lower level of the controllers is discussed for the

leader and follower satellites pair to maintain and track the desired relative position and

attitude. First, in this Chapter, controllers are designed for orbit maintenance; then in

Chapter 4, attitude controllers and coupled attitude and orbit controllers are designed.

In the literature, numerous controllers have been proposed for the relative position

and attitude for the leader-follower formation flying system [21, 22, 26, 31, 40, 41, 47,

49, 58, 64, 65, 66]. Roughly, these controllers can be classified into either exact-model-

based or uncertain-model-based controllers.

In the exact-model-based controllers, Linear Quadratic (LQ) and Proportional-

Integration-Differential (PID) are often used. It appears that the LQ controller is the first

controller to have been applied to the formation flying system [21, 22, 26, 40, 58]. Based

on the Clohessy-Wiltshire (C-W) model, Vassar and Adams [21, 22] proposed the









continuous LQ controller in order to keep the relative position within the formation

system. In order to design a controller for this relative position model that uses a pulse

based thruster, Kapila [26] made some variations to the C-W model. One was discretizing

the continuous model, and the other is using Av instead of normalized forces. Through

these changes, the C-W equation was developed as xl, = ODx + FAv, where i and i +1

represent discrete time steps. D and F are constant matrices. Based on this discrete

model, the pulse based controller was designed.

Other exact-model-based controllers found in the literature are the PID class of

controllers. For example, based on the C-W model, Robertson [40, 58] proposed a PD

controller augmented with a thruster mapper. The controller has a format as shown in Eq.

(0.1), where g(y) is the thruster mapper and is defined as shown in Eq. (0.2). Kp and Kd

are the proportional and differential gain and k is the state error. y is a performance

criterion value and fh,,,,,, is a constant thruster force level. Additionally, Robertson [58]

proposed a bang-off-bang controller, which is also based on the C-W model and has a

similar idea to the one presented in reference [40].

u= g(y)(K + Kd,) (3.1)

+2 y<1.5fthuster
+1 y > 0.75fthrust
g(y)= -1 y <-0.75fthrr (3.2)
-2 y < -0.15fhrt
0 oh'I/I Il it'

All of the above exact-model-based controllers are based on the linear C-W

model, which has some limitations in representing the true motion. They all assumed that

the system model is well defined and all the states are known. Since, this assumption is









not always realizable, and the controllers based on this assumption are sensitive to the

spacecraft's parameter uncertainties, unexpected disturbance, and the initial

attitude/position rates, the resulting performance will be less then desirable. Therefore,

uncertain-model-based controllers [31, 41, 47, 49, 64] have appeared in the literature

recently. For example, adaptive, Lyapunov-type and sliding mode controllers are now

widely used for the system with model uncertainty.

Bernstein [47] proposed an adaptive asymptotic tracking strategy for the

spacecraft attitude of the formation system in the case that the spacecraft's inertia matrix

is not known. By re-writing the attitude governing equation, the large-angle-tracking

problem becomes an attitude-keeping problem near the desired point. Using this proposed

adaptive law, the tracking error can be controlled to zero. What's new in Bernstein's

contribution is that, at the same time, the spacecraft's inertia matrix can be predicted.

However, in Bernstein's work the attitude and position are assumed to be uncoupled,

whereas in this Chapter, the controller will be designed for the coupled system.

A Lyapunov-type controller is another typical method used for the uncertainties

model, which is mainly applied to the relative position control problem. For instance,

based on the circular reference orbit, Kapila [31] modeled the uncertainty as

4 =[mF,mFMeG,mF /mL,,ad,axdadT where m and a represent the mass and the

normalized acceleration of a satellite separately. Me and G are the Earth mass and the

gravitational constant. Subscripts L and F are the leader and follower, whereas, x, y,

z and d represent the coordinates and disturbance. Then the model error is = 4-4,


where 4 is the predicted model parameters. Accommodating this uncertainty, the









Lyapunov-type controller is designed. This controller guarantees global asymptotic

convergence of the position tracking error in the presence of unknown, slow varying

satellite masses, disturbance forces, and gravity forces. Nevertheless, this method need to

be extended for the nonlinear coupled system.

Kim and Choi et al [49, 64] used the sliding model control method in their work

on spacecraft attitude and the smart structure control. Using this sliding mode controller,

the uncertainties related to the natural frequency and external disturbances of the

materials can be rejected. The sliding model will be demonstrated in this Chapter to have

the capability of controlling the formation system well.

All of the above uncertainty-model-based controller examples demonstrated the

capability of system control under parametric uncertainties, such as mass, inertia, and

gravitational constant. However, there are some other uncertainties and limitations that

have not been considered; for example, un-modeled dynamics, actuator saturation and

actuator resolution. Furthermore, the full nonlinear model should be used instead of the

linear C-W equation in order to obtain the precision formation.

For the formation flight examples considered above, the dynamics models are

developed based on the leader-follower (or hierarchical) structures. There are, however,

some other ways to address formation systems. For example, Hadaegh and Goulet [65,

66] considered the thruster saturation and treated the formation system as a non-

hierarchical structure. An adaptive saturated controller is designed based on this

assumption. This adaptive controller continuously updates the estimated satellite mass

and the mass is regarded as the uncertainty for this model.









In this Chapter, the controllers based on the generalized model derived in Chapter

2 are designed. First, a performance criterion is discussed in Section 3.2. Then, in Section

3.3, in order to show the advantages of the sliding mode controller based on the full

nonlinear model, controllers based on the following methodologies are develop: LQ

method, /u synthesis, feedback linearization, and sliding mode controllers. Next, a

comparison of the performance of these controllers is conducted in Section 3.4. The

chapter concluded with a summary of the results obtained for the orbit control problem.

3.2 Performance Criterion Of The Leader-Follower Pair

In this Section, the control performance criterion is given. In this Chapter, only

the lower level leader-follower pair's control problem is addressed. Therefore, the

criterion of the controller for the formation flying system can be simplified to the

criterion of the controller for the leader-follower satellites pair. In most of the small

satellite formation flying missions, such as UoSAT-12 and ACE, the relative position

needs to be kept within a certain range and the attitude of the satellites needs to be

tracked or maintained cooperatively. Therefore, the criterion of the controller is to

minimize the above two branch errors (relative position and relative attitude errors) with

minimum control efforts.

Performance Criterion: The actual and desired relative position between the

leader and follower satellites are p and p Then the error of the relative
.ACTUAL -DESIRED

position is p =p p DESRED The error of the attitude is modeled as
S[EL,AZ] in the small Euler angle sense. Let the generated control, i.e., force

ERROR = [EL,AZ]T in the small Euler angle sense. Let the generated control, i.e., force









and torque commands, to be u. Therefore, the criterion of the formation system is shown

as Eq. (3.3), where Q Q and R are weighted matrices with proper dimensions.
=1 =2

J= l j oQ p +OERRoRQ O o+uT Ru)dt (3.3)
RROR =1 =2 --
0 ERROR

In the following sections, linear and nonlinear controllers are developed based on

minimization of the above criterion. Detailed performances for each method are

explained in the corresponding section.

3.3 Relative Position Controller Design

In this section, four control strategies are developed and implemented; they are

LQR, u/, feedback linearization, and sliding mode control. First, based on the linear

dynamics model (C-W equation), the LQR and /u controllers are designed with

uncertainties. Then, the feedback linearization method together with the /u synthesis is

studied for the nonlinear model. Finally, the nonlinear control strategy, sliding mode

control, is designed. All the controllers are assumed to have the full state feedback. The

simulation results, discussion, and comparison are given in Section 3.4.

3.3.1 Linear-Quadratic State Feedback Regulator (LQR)

The LQR method is used to design the optimal controller for the linear system in

H2 sense.

3.3.1.1 LQR

Theory 3.1 ([67], p273): Assume the linear state space system is modeled as Eq.

(3.4), where A, B, C, and D are constant matrices. The performance criterion is shown

in Eq. (3.5) where Q > 0 and R > 0 are the weighted matrices with proper dimensions.









x= Ax+Bu
-- (3.4)
y = Cx+ Du


J = x Qx+u Ru)dt (3.5)
2 -

If matrix pairs (A, B) is controllable and (A, C) is observable, then there is one


and only one optimal controller as shown in Eq.(3.6), in which P satisfies the Riccati

equation (3.7).

= -R BT Px(t) (3.6)


PA+AP- PBR B TP+Q=0 (3.7)


3.3.1.2 LQR synthesis model

For the linear circular reference orbit, the relative position dynamics model can be


written as the C-W equation as shown in Eq. (3.8), where a, /, and n = fL are
\a

the semi-major, gravitational parameter of the central body (Earth), and mean motion of

the leader's orbit respectively. The state vector includes the relative position and velocity

expressed in LVLH as xn =[x,y,z, x, y, z]T, where subscript n denotes the nominal

1
model. The actuator, which is assumed to be is not easy to be synthesized in the
s+l

LQR controller's design, therefore, the system equation is rewritten to include the

actuator model before the design. First, the actuator's transfer function is written in the

system equation form as in Eq. (3.9), where the subscript a represents the actuator. After

synthesis of the actuator system equation with the C-W model as shown in Figure 3.1, the










combined linear system equation is as Eqs. (3.10) and (3.11), where the state vector is

x=[x,x, '.


0 0 0
0 0 0
0 0 0
3n2 0 0
0 0 0
0 0 -n2
0 0 -n2


1 0 0
0 1 0
0 0 1
0 2n 0
-2n 0 0
0 0 0


Xn +


0 0 0
000
0 0 0
0 0 0
1 0 0
0 1 0
0 0 1
001


a = -x +u

ya = xa



- Actuator C-Wmodel


Figure 3. 1 Combined state space



0 0 0 1 0 0 0 0
0 0 0 0 1 0 0 0
0 0 0 0 0 1 0 0
3n2 0 0 0 2n 0 1 0 (
0 0 0 -2n 0 0 0 1
0 0 -n2 0 0 0 0 0
0 0 0 0 0 0 -1 0 (
0 0 0 0 0 0 0 -1
0 0 0 0 0 0 0 0 -
Ax+Bu


(3.8)


(3.9)


0 0 0
0 0 0
0 0 0
0 0 0
0 0 0
0 0 0
1 0 0
0 1 0
0 0 1
000
000
000





100
010
001


u
(3.10)









1 0 0 0 0 0 0 0 0 000
0 1 0 0 0 0 0 0 0 0 0 0
0 0 1 0 0 0 0 0 0 0 0 0
y-- _x+ u
S0 0 0 1 0 0 0 0 0 0 (3.11)
0 0 0 0 1 0 0 0 0 0 0 0
0 0 0 0 0 1 0 0 0 0 0 0
= Cx+Du

The controller's closed loop objective is to track and maintain the desired states.

As such, r = [r,r, r, r, r, r, r6]T is assumed to be the commands or desired states, and the

error e = r- n _- is aimed to zero. In order to design the controller based on LQR

synthesis, a new state is chosen as in Eq. (3.12), and new system equation is derived as

Eq. (3.13).

S= [x, edt]T (3.12)


A 09 9 B 0 9x6
x= x + u+ r (3.13)
-C 06x6 06 63 6 6

In order to minimize the errors that are last three elements of the state, weight

functions for the state and control are chosen by try and error as in Eqs. (3.14) and (3.15).

Using MATLAB LQR command, the feedback gain for the system can be obtained as in

Eq. (3.16), where K represents the column vectors from columns i to j. K K A
=(,']J) =J =A

and K are the gains for nominal model, actuator model, and integrator model
INT

respectively.

0.035I 1 0.1
S 03512x12 12x3 (3.14)
0312 I33










F 1007
R= 0 1 0 (3.15)
0 0 1


K=[K K ,K
-A NT (3.16)
=[ (,16)' (,79)'= (,1015)

3.3.1.3 LQR analysis model

Figure 3.2 gives the analysis model of the LQR controller. The input is the desired

relative position (i.e., commands), and the outputs are the relative position and velocity in

LVLH coordinates and actuator forces. P and Av are the nominal C-W and actuator

models. The integral model 1 is used to reduce the steady state error. Furthermore, the
S

uncertainty model is shown in Figure 3.3. P represents the C-W model with a

uncertainty in the leader orbit's mean motion that has 200% uncertainty. The un-modeled

dynamics related to the nominal model and actuator are bounded by A, and A2, where

50% and 20% uncertainties are assumed.



























Figure 3. 2 Analysis model of LQR controller without uncertainty


Figure 3. 3 Analysis model of LQR controller with uncertainty


























Figure 3. 4 Typical loop gain and phase margin model



In the classical control design, the gain and phase margins are used to study the

robustness. The gain margin is the amount of gain increase required to make the loop

gain unity at the frequency where the phase angle is -1800. The phase margin is the

difference between the phase of the response and -1800 when the loop gain is 0 dB. It is

generally found that gain margins of six or more combined with phase margins of 450 or

more result in reasonable trade-offs between bandwidth and stability [46]. Each loop in

the analysis model needs to be tested, and for simplicity, only one loop gain and phase

margin's model diagram is shown as in Figure 3.4. In this model, the loop gain and phase

margins from command the r, to the output, relative position x, is tested. The control

gains, time response and robust analysis based on LQR will be shown in Section 3.4.

3.3.2 H //u Controller

3.3.2.1 Structured singular value (u/) theory

Hl / / is a tool used to design robust controllers and to analyze the robustness

and performance of a system subjected to uncertainties, especially the structured









uncertainties. The structured uncertainty represents the normalized uncertainty, which has

a fixed structure with scale parts and matrix parts as shown in Eq. (3.17). 5, is a scale, In

is a r x r identity matrix, and A, is a square matrix. The H. / / controller searches only

on the structured uncertainty, so it is not as conservative as based on the small gain

theory. /u is defined as the inverse of the largest uncertainty, which causes the system to

be unstable as shown in Eq. (3.18), where a(A) is the singular value of the uncertainty

A, and M is the open loop model. If the /u value of the closed loop system is less than

1, the controller is robust with respect to the structured uncertainties.

A = {A = diag[8,I,,, 8,12,.,1 ,, .. ,aF ]:
5, GeC,n E C"J m 3 <1, A <1, (3.17)

i [1,s], je[1, F i|A|
1
P -= (3.18)
min{(a ): det(I-MA) = 0}
GA

/u synthesis is very difficult to be implemented, and D-K iteration is normally

used [68, 69]. The basic ideas of the D-K iteration are first scaling the open loop plant by

a state-space D, then compute a optimal H. control gain K for the scaled plant.

Second, compute set of optimal scaling for the closed-loop system under the fixed K to

reduce the /u bounds. Iterating these two steps until the peak /u value does not change

much. The detailed information of / synthesis theory and D-K iteration can be found in

Lind [68] and Matlab [69].

3.3.2.2 /1 synthesis model

First, the parametric uncertainty model, which has 200% uncertainty in the mean

motion of the leader, is derived. Two correlated parametric uncertainties are considered.









One is the mean motion of the leader c2 20 (1+ W2A2), where c2 = n. The other is the

square value of the mean motion c, = c,1(1 + WA1), where c0 = n2 c20. A, and A2 are

normalized uncertainties. W = W22 = 4 are uncertainties' magnitude in this dissertation.

Substituting these two parametric uncertainties into Eq. (3.8), the system equations can

be rewritten as Eqs. (3.19) through (3.21). Therefore, the open loop parametric

uncertainty state space model can be written as Eq. (3.22), which has the input-output

relations as shown in Figure 3.5. wo, i e [1, 4] are the exogenous signals, which includes

disturbance and noise, u,,ie[1,3] are the control inputs, z,,ie[1,4] are the outputs

related to o,, and output y, i e [1,6] are composed of the relative position and velocity.

X4 = 3cox, + 2c2X5 +u +3c1oWAlxl +2c20W2A X5
=3cox, + 2c20X5 +u1 +ZA1, +Z2A2 (3.19)
= 3cox, + 2c20X5 +U1 +01 +02

x5 = -2c2o4 + u2 2c20W1x4
= -2c20X4 +2 + Z3A1 (3.20)
= -2c20x4 + u2 +)3

6 =-Co10X+3 3 C10W2A2X3
= -CoX3 + U + Z4A2 (3.21)
= -CIX3 + 3 +04















3c10
0
0
3c10Ff
0
0
0
1
0
0
0
0
0

A
=6x6

_10x6


0
0
0
0
0
0
0
0
0
0 -
0
1
0
0
0
0

B
6x7
o

























=67


0
0
0
0
0
-C10
0
0
0


0
0
1
0
0
0

1


1
0
0
0
-2c20
0
0
0
-2c20W
0
0
0
0
1
0
0
o






















0


0
1
0

2c20
0
0
0

2c20W2
0
0
0





1
0


(3.22)


D
=10x7


Z1
Z2
Z'
z,
Yl
Ya

Y4
Y5
Y6


Figure 3. 5 Open loop model with parametric uncertainty


The /j synthesis model is shown in Figure 3.6, where A= diag{A,A 2, ,A2} is


the parametric uncertainty. A and A are normalized un-modeled uncertainties
associated with the C-W equation and the actuator, respectively. W and are the
associated with the C-W equation and the actuator, respectively. Wm and WA are the


"1,

ats
09t.

u,
EdP
";
"i










uncertainties magnitude with values of 50% and 20%. The measurement noise is assumed

to be W = 0.021 The weighting function for the performance and control efforts are
-Pn =6x6

O.1s 0.1(s+1)
selected as W =0 I and W =. (s I respectively. The actuator is
=perf s+0.1=66 =u s+10 =3x3

1
assumed to be the same as in LQR synthesis, which is -- Figure 3.7 gives the input-
s+l

output relations of this synthesis model. The differences between the synthesis model and

the open-loop model are commands r, noise n, performance e, and control effort u.

The feedback gain K between u and y are found through D-K iteration commands in

MATLAB.


Figure 3. 6 Synthesis Model











01 Z,




r(6} Q e_,(6}

I(6M e_ (3)
u(3) ypN)






Figure 3. 7 Synthesis model input-output relation



3.3.2.3 u/ analysis model

-xy,z}
Ti( n6)i w
r 6! (Xyz}


{xy,z,



Figure 3. 8 Analysis model for nominal Figure 3. 9 Analysis model for actual
system system



After using the synthesis model to design the controller, the analysis model is

constructed to evaluate the controller for the time response, frequency response, and

robustness analysis, respectively. The analysis model for the nominal system and actual

system are shown in Figure 3.8 and Figure 3.9. These two diagrams are similar to the

ones in LQR analysis model and not discussed here.

3.3.2.4 NP, RS, and RP

In H. /I/ type controller design, NP, RS, and RP are used to study the robustness


and performance. RS: robust stability, family of models with uncertainty, which is stable.









NP: nominal performance, family of models without uncertainty, which achieves the

desired goals. RP: robust performance, family of models with uncertainty, which

achieves desired goals [68]. After obtaining the feedback gain K, the closed-loop model

with uncertainty is formed as shown in Figure 3.10. Furthermore, the equivalent models

S S
can be obtained as in Figure 3.11 and Figure 3.12. Here, S = 2 = FL (G, K) is
S S
=21 =22

the lower Linear Fractional Transformations (LFT) of the models G and K with proper

dimensions corresponding to the input and output.






G

y(6) I I I 1 33) 6 S
e-=[e,{(6);{3}] d=[r_{6};6}



Figure 3. 11 Equivalent model 1
Figure 3. 10 Closed-loop model with
uncertainty






Z, 6 Z1 1 d=r 6)13



e = [e{ e.(3] I d=[L(6);n(6)


Figure 3. 12 Equivalent model 2