UFDC Home  myUFDC Home  Help 



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 FitzCoy, 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 goodspirited 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 2Body 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 (ClohessyWiltshire 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 LeaderFollower Pair......................69 3.3 Relative Position Controller Design ................ ... ........................70 3.3.1 LinearQuadratic 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 CW Model................... ............ 96 3.4.4 Time Response Based on CW 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 21 Different cases been simulated............. ......... .................... 29 22 Classical elements for an Earth satellite........... ............ ................59 31 G ain and phase m argins............... ................ .............. ....... ......95 32 Robustness performance and controller order comparison.......................96 33 Fuel consum ption................................... ..................... .......... 98 34 Fuel consum ption............................... ................ ..............101 4 1 C ases ................... ................... ..... .................. .... ....... 107 42 C ases ................... ................... ..... ................... .... ....... .110 43 C ases ............. ............................. .................... ... ......... 113 44 GA parameters' definition and the data used..................................... 132 45 Simulation results of GA + SMC + BLC............... .......... ..............133 51 The initial orbit information of the formation satellites..................... 139 52 Orbit information of R1 and R2 in figures 3 and 4 ........................... 142 53 Formation maintenance controllers for follower 1 and follower 2.............147 54 Formation reconfiguration controllers for redundant 1 and redundant 2........147 55 Configuration code for inplane 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 microsatellite 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 inputoutput relation ........................................ ........................ 81 3.8 Analysis m odel for nom inal system ....................................................... 81 3.9 Analysis m odel for actual system ............................................................................81 3.10 Closedloop model with uncertainty..................... .. ................ ................. 82 3.11 Equivalent model 1 ................................. .. .. ........................ 82 3.12 Equivalent m odel 2 ................................... ............. .. ............82 3.13 Closedloop analysis model based on the feedback linearization.................................89 3.14 Closedloop 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 xdirection. ................................. ............................... 134 4.59 Control force in xdirection. ............................................... .............................. 134 4.60 The relative position error in ydirection. ................................. ............................... 134 4.61 Control force in ydirection. ............................................... .............................. 134 4.62 The relative position error in zdirection. ........................................ ............... 134 4.63 C control force in zdirection ......... ................. ................... .................. ............... 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 Leaderfollow er Structure................................................................. ............... 137 5.3 A b solute and interconnection ........................................... ........................................ 137 5 .4 Su b form action ............................. ........................................................... ............... 13 7 5.5 3satellite 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 550010000 se c o n d s ...................................... ............................... ................ 1 5 0 5.12 Relative position differences in z local frame, (a) 0800 seconds, and (b) 56506500 seconds ................................................................ ... ..... ......... 151 5.13 Force commands in x local frame. (a) whole range, and (b) 560010000 seconds ......152 5.14 Force commands in y local frame. (a) whole range, and (b) 56009000 seconds ........152 5.15 Force commands in z local frame. (a) 0800 seconds, and (b) 59006300 seconds .....153 5.16 Elevation error angle. (a) whole range, and (b) 567015000 seconds ..........................154 5.17 Azimuth error angle. (a) whole range, and (b) 59007000 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 FitzCoy Major Department: Mechanical and Aerospace Engineering Due to the potential advantages of the microsatellites 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 multisatellites 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 ClohesseyWiltshire 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 attitudeorbit relative motion and to design controllers based on this model. The developed generalized dynamics model is based on the leaderfollower 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 twotiered 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 + Selfcontrolling vehicles eliminate the need for extensive ground support With respect to the conventional single large satellites, a formation flying system may provide numerous advantages [58]. 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 lowEarth 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 microsat will be about 100 Kg, and that of a nanosat will be approximately 20 Kg. Therefore, at least 50 nanosatellites 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 coobserving programs can be conducted without using extensive ground support. For instance, in the leaderfollower 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 Xray 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 groundtomultiplesatellite communications path is more complex than the groundtosinglesatellite 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, 1012], 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, lowmass 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, EO1 is actually the first prototype in the New Millennium Program of a formation flying system. EO1 works with the LandSat7 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, EO1 and LandSat7 will image the same point in order to get the piecewise continuous data. The EO1 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 formationflying concept because EO1 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 EOl's small satellite. The mass of the small satellite is 90 Kg. Figure 1. 3 The satellite for EO1 The objective of the EO3, which is also called GIFTIOMI [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 20052006. The mass of the instrument is about 100 Kg. Satellite Technology 5 (ST5) 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 microsatellite and the communication architecture. The righthand 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 Xray, ACE [15], UoSAT12, 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 subsystems. The ORION microsatellite 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 ClohesseyWiltshire (CW) equations as the dynamic model upon which the controller is designed [1930]. CW 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, noncircular, 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 leaderfollower 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, 3234]. 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 multitiered control architecture that has been used in the smart highway system [35]. In Stadter [36], the highlevel 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 stateoftheart of the formationflying 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 leaderfollower 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 sixorbital elements will be discussed. Chapters 3, 4, and 5 address the control strategies for the formation flying system. A twotiered 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 reconfiguration 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 ClohesseyWiltshire (CW) equation [1924, 2630, 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 inplane and outofplane 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 CW 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 CW 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. Lyapunovtyped 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 Tperiodic 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 CW 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 openloop simulation environments for spacecraftrelated 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 Nbody problem is simplified to the 2body 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 Nbody 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 nongeostationary), 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 intersatellite 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 microsatellites 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 intersatellite attractions are even smaller than those nonspherical 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 Earthsatellite 2body relative motion problem as shown in Eq. (2.5). Here 9ad includes the normalized control efforts, disturbance and Nbody attractions. R + =R = ad (2.5) 2.3 General 2Body Relative Motion Models In this section, the general 2body 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 2body 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 EarthCenteredInertial (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 righthand 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 2body 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 unmodeled 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 P2 _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 22 + 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 timedependant solutions of the classical 2body 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 2body 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(1e) /a(1e )(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 p2 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(1e2), 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 semimajor 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) Rx 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) Rn 2 cos Rv 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 cos1 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 Re 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 nonzero 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 21 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 2O ,,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 21. 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 2norm difference percentage shown in the results (Figure 2.13) is on the order Table 21 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 S50 620 60 x axs in Local Frame (kim) Figure 2. 1 Relative position in 3D x 1016 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 21 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 CW equation and the inplane nonlinear circular reference orbit equation. The third one is the inplane nonlinear elliptic reference orbit equation. Here inplane represents the case where the leader and follower satellites are in the same plane. 2.3.3.1 Linear model (ClohessyWiltshire equation) If the reference orbit is a circular orbit, the general 2body relative motion dynamic model derived in Section 2.3.1 can be simplified to a linear system, which is the wellknown ClohessyWiltshire (CW) 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 2body 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 CW 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 inplane and the outofplane motion, nonlinear coplanar 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 inplane 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 inplane elliptic reference orbit relative dt df dt' motion can be rewritten 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 inplane 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(1e ) (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 CW equation is are shown by simulation. 2.3.4 Simulation Results and Discussion (Point Mass Case) Simulation results show the difference between CW equations and the generalized relative position model. Results demonstrate the necessary of using nonlinear model instead of the CW linear model in the formation flight system. The inplane circular and elliptic reference orbit model are not compared, since it is identical to the generalized relative position model under the inplane 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 CW equation and the generalized relative dynamics model. In Figure 2.18, the difference between two models is shown. The difference is calculated by the 2norm between these two models at each time second. CW 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 CW 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 S20 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 EO1 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 realtime 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 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 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 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 123 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 nonzero 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 semimajor 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 ST5 [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 0degree 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 nonuniform 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 2norm 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(104) 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 14 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 22. 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 semimajor 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). acosEe R=C bsinE (2.93) 0 a f sinE/(1 e cosE) = C bfcosE(1ecosE) (2.94) 0 E le 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 2body 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(1e) 1, 80 0 OR 1 8= a(1 e) m2 1n2 OR Sa na m 8e nL 2 S=fb/(le) 12 R1 = a(1 e) sin o m, 8i n,/ LR i 3bnt 12 S=(1e) mn,  m Oa 2q In2 = b /(1 e) m, LiJ 1 2 1 12 /F 8l 3f t /b = j bcos0/(1e) m3 3 bf m 2 (2.107) Oi 3 a 2(1e)2 1 2a( e) 2 L/3 1, /2 L['2 1' 8e b(1e) 2 A (1e)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  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 reconfigure the formation. The controllers for the formation system can be treated as a twolevel architecture, in which the higher level dealing with the systematic control and the lower level related to the individual leaderfollower 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 leaderfollower formation flying system [21, 22, 26, 31, 40, 41, 47, 49, 58, 64, 65, 66]. Roughly, these controllers can be classified into either exactmodel based or uncertainmodelbased controllers. In the exactmodelbased controllers, Linear Quadratic (LQ) and Proportional IntegrationDifferential (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 ClohessyWiltshire (CW) 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 CW model. One was discretizing the continuous model, and the other is using Av instead of normalized forces. Through these changes, the CW 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 exactmodelbased controllers found in the literature are the PID class of controllers. For example, based on the CW 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 bangoffbang controller, which is also based on the CW 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 exactmodelbased controllers are based on the linear CW 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, uncertainmodelbased controllers [31, 41, 47, 49, 64] have appeared in the literature recently. For example, adaptive, Lyapunovtype 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 rewriting the attitude governing equation, the largeangletracking problem becomes an attitudekeeping 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 Lyapunovtype 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 = 44, where 4 is the predicted model parameters. Accommodating this uncertainty, the Lyapunovtype 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 uncertaintymodelbased 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, unmodeled dynamics, actuator saturation and actuator resolution. Furthermore, the full nonlinear model should be used instead of the linear CW equation in order to obtain the precision formation. For the formation flight examples considered above, the dynamics models are developed based on the leaderfollower (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 LeaderFollower Pair In this Section, the control performance criterion is given. In this Chapter, only the lower level leaderfollower 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 leaderfollower satellites pair. In most of the small satellite formation flying missions, such as UoSAT12 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 (CW 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 LinearQuadratic 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 CW equation as shown in Eq. (3.8), where a, /, and n = fL are \a the semimajor, 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 CW 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 CWmodel 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 CW 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 CW model with a uncertainty in the leader orbit's mean motion that has 200% uncertainty. The unmodeled 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 tradeoffs 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 iA 1 P = (3.18) min{(a ): det(IMA) = 0} GA /u synthesis is very difficult to be implemented, and DK iteration is normally used [68, 69]. The basic ideas of the DK iteration are first scaling the open loop plant by a statespace D, then compute a optimal H. control gain K for the scaled plant. Second, compute set of optimal scaling for the closedloop 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 DK 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 inputoutput 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 unmodeled uncertainties associated with the CW equation and the actuator, respectively. W and are the associated with the CW 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 openloop model are commands r, noise n, performance e, and control effort u. The feedback gain K between u and y are found through DK 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 inputoutput 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 closedloop 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 Closedloop 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 