UFDC Home  myUFDC Home  Help 



Full Text  
ARTIFICIAL POTENTIAL FUNCTION GUIDANCE FOR AUTONOMOUS IN SPACE OPERATIONS By ANDREW R TATSCH 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 2006 Copyright 2006 by Andrew R Tatsch Dedicated to My parents who brought me into this world and My friends, family, and colleagues who opened my eyes to it. ACKNOWLEDGMENTS To my advisor, Dr. Norman FitzCoy, I would like to convey my gratitude. Few graduate students have the opportunity to work with, not for, an advisor that encourages freedom and creativity in the search for a pragmatic solution and uniquely challenges them beyond what they think is possible. He also taught me many lessons beyond the curriculum and for those I am grateful. Without his guidance this dissertation would never have come to fruition. Many thanks also go out to my supervisory committee, Drs. Rick Lind, Andrew Kurdila, Gloria Wiens, BJ Fregly, and Stanley Dermott, for their invaluable assistance and patience. Lastly, I would like to acknowledge my brother J Bradley Tatsch for everything that he has done to make this dream a reality. Without his assistance my undergraduate and graduate experiences would have never gotten off the ground. I am forever grateful for all that he has done for me. TABLE OF CONTENTS page A C K N O W L E D G M E N T S ................................................................................................. iv LIST OF TABLES ............ ........... .. ......... ................. ........... vii LIST OF FIGU RE S ............... .... .......... ............... ...... ......... .. .......... viii A B ST R A C T ................. .......................................................................................... x CHAPTER 1 IN TR OD U CTION ............................................... .. ......................... .. State of the Art for InSpace Operations (ISO) Functionalities .............. ...............4 R research Scope ...................................... ............................. ................. 11 D issertation O outline .................. ............................................. .............. .. 12 2 TRAJECTORY DETERMINATION: ORBIT MECHANICS APPROACH........... 14 M idcourse Trajectory D eterm nation .............................................. .....................14 Lam bert's Problem ............................................... ........ ............... 16 M maneuver Command Determination .... .......... ........................................ 19 N um erical Exam ple: M idcourse................................... .................................... 21 Endgam e Trajectory D eterm nation ...................................... ........... .................. 22 CW Maneuver Command Determination............ ...... ..............24 N um erical Exam ple: Endgam e.................................... ..................................... 25 Alternative Orbit Mechanics Based Methods..........................................................26 C o n clu sio n s.................................................... .................. 2 8 3 ARTIFICIAL POTENTIAL FUNCTION POSITION GUIDANCE.........................29 R ob ot M option P planning .................................................................... .....................2 9 Artificial Potential Function Position Guidance...................................................32 Example: Spacecraft Unconstrained Rendezvous ............................................... 35 O b struction s .................................................................................. ......... 37 Numerical Example: PathConstrained Rendezvous..............................................42 C o n clu sio n s..................................................... ................ 4 4 4 ARTIFICIAL POTENTIAL FUNCTION GUIDANCE: A NEW APPROACH.......45 Dynamic Artificial Potential Function Obstacle Avoidance Algorithm...................45 A attractive Potential .................. .............................. .... .. .. ......... .... 46 R epulsive Potential .................. ............................. .. ... ... .. ........ .... 47 T otal F orce Interactions........................................................................... 49 Sim ulation R results ............................................... ........ .. ............ 50 C constant velocity scenario ........................................ ....................... 50 O norbit rendezvous scenario.................................... ........ ............... 52 Optimized Harmonic Attractive Potential ...................................... ............... 54 C o n c lu sio n s..................................................... ................ 5 7 5 ARTIFICIAL POTENTIAL FUNCTION ATTITUDE GUIDANCE.....................60 A attitude E quations of M otion ......................................................... .....................63 Continuous APF Attitude Guidance .................................. ............... ............... 64 Stab ility A n aly sis............ ...... ............ ............................................ .... .... ..... 6 7 Stability of Dynam ical System ................................ ......................... ........ 72 O ptim ization ...................................... ............................ ............... 74 N um erical Exam ple .............................................. .... .... .. ........ .... 75 Im pulsive A PF A attitude G uidance................................................... ......... .......... .. 79 Numerical Example: Impulsebased APF Attitude Guidance.............. ................ 82 C o n c lu sio n s.................................................... .................. 8 5 6 CONCLUSION AND FUTURE WORK ........................................ .....................87 C o n clu sio n ...................................... ................................. ................ 8 7 F u tu re W o rk .......................................................................................................... 8 8 L IST O F R E FE R E N C E S ....................................................................... ... ................... 90 BIOGRAPH ICAL SKETCH ...................................................... 96 LIST OF TABLES Table page 1 Orbit parameters for rendezvous simulation. ................................ .................52 2 G enetic algorithm statistics ............................................. ............................. 57 3 Mean and standard deviations for the two parameters in the GA optimization .......57 4 G genetic algorithm statistics ........................................................... .....................74 5 Mean and standard deviation for the six parameters in the GA optimization ..........75 LIST OF FIGURES Figure page 1 M odel of the ROTEX m anipulator................................ ........................ ......... 5 2 Artist's rendition of ETSVII components ..... ......... ....................................... 5 3 CM U 's Skyw orker prototype........................................................... ............... 6 4 N A SA 's Robonaut ....................... ................................................. .6 5 University of Maryland's Ranger ............ ........ ...........................6 6 AERCam Sprint.................. ......................... ........... .... ............. 7 M ini A ER Cam .................................................... 9 8 Artist's depiction of the Orbiter Boom Sensor System .................. .................. 9 The SCAMP in the neutral buoyancy tank..............................................9 10 Schematic of Lambert's Problem.................... ................. ... ............... 15 11 Transfer trajectory from midcourse trajectory determination aglorithm................23 12 D epiction of orbital rendezvous ........................................ .......................... 23 13 Endgame maneuver trajectory in the xyplane (Rbar Vbar plane).....................26 14 Glideslope solution for 5 impulse rendezvous. .............................. ................27 15 Planar rendezvous in the CW frame under APF guidance............................. 36 16 Control efforts in x and y directions for the rendezvous........................................36 17 Difference between IC potential and potential on surface of obstructions ..............40 18 Contour plot of the potential function.......................................................... 40 19 Zoom ed view of contour plot......................................................... ... ............ 41 20 Stem plot of obstruction surface potential error................................................41 21 Trajectory traversed by robot for the pathconstrained rendezvous example ..........42 22 Control action in the xdirection for the constrained rendezvous ..........................43 23 Control action in the ydirection for the constrained rendezvous ..........................43 24 2D Constant velocity simulation results........................................ ............... 51 25 Planar onorbit rendezvous results. ........................................ ....... ............... 53 26 Norm of applied chaser accelerations. ........................................ ............... 54 27 Optimized sinkvortex pair vs. CW based vector field.............. .............. 58 28 Inequality of E q. (64) ...................... .................. ................... .. ...... 72 2 9 B eh av ior of q4 ................................................. ............... ....................73 30 Time response of vector part of quaternion .........................................................76 31 Tim e response of scalar part of quaternion................................... .................77 32 Tim e response of angular velocity ........................................ ....................... 78 33 Tim e response of control effort.................................................................... ..... 79 34 Components of the vector part of the quaternion vs. time .................. .............82 35 Components of the angular velocity vs. time .......................................................83 36 Scalar part of the quaternion vs. tim e.................................................................... 83 37 Angular velocity impulses vs. time for the numerical example..............................84 38 Scalar part of the quaternion vs. time for minus one convergence case ..................85 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 ARTIFICIAL POTENTIAL FUNCTION GUIDANCE FOR AUTONOMOUS IN SPACE OPERATIONS By Andrew R Tatsch May 2006 Chair: Norman FitzCoy Major Department: Mechanical and Aerospace Engineering For most inspace operations (ISO), e.g., assembly, repair, refueling, the target location and its orbit parameters are known a priori. For instance in the case of an assembly mission, the general location of the construction site and the parts depot are known parameters in mission planning. Utilizing this known information, we propose a novel approach to collisionfree operations for a heterogeneous team of spacecraft robots performing the ISO mission based on artificial potential function (APF) guidance. In order to utilize APF guidance two prominent deficiencies of APF methodologies need to be circumvented: premature termination due to local minima in the operational space and suboptimal performance. It is widely known that utilization of harmonic potential functions eliminates the introduction of local minima in the operational space; thus we employ harmonic functions to conquer the first of the aforementioned deficiencies. In order to surmount the suboptimal performance, we propose the use of a sum of harmonic function primitives (simple functions such as sinks and vortices in the hydrodynamic analog) to shape the velocity field of the operational space to harness the underlying dynamics of the system (i.e., orbit mechanics). Furthermore, given the location and orbit parameters of the target site, a random heuristic search algorithm can be employed to optimize the velocity field shaping for maximizing the performance gains achieved from harnessing these underlying dynamics. CHAPTER 1 INTRODUCTION I believe that this nation should commit itself to achieving the goal, before this decade is out, of landing a man on the Moon and returning him safely to the Earth.1 President John F. Kennedy, May 25, 1961 Just as the speech given fortythree years ago on a Thursday in late May by President John F. Kennedy presented the National Aeronautics and Space Administration (NASA) a new directive, on January 14, 2004, President George W. Bush delivered a speech titled "Vision for Space Exploration,"2 whose intent was to again unify NASA around a set of common goals. These goals are centered around sending man and/or robotic missions to the Moon, Mars, and deeper into our solar system. Implicit in the robotic missions proposed by this new space exploration mandate are the utilization and development of autonomous space robotic systems that "will serve as trailblazers, advance guard to the unknown."2 Based on these new directives, autonomous space robotic systems will be imperative for the success of NASA's new vision. Space robotics, in the context of this work, are robotic systems whose applications can be roughly dichotomized into two broad classes of missions: inspace operations (ISO), and planetary surface exploration (PSE). Both of these mission profiles share common functional requirements and possess functional requirements unique to each mission. The functionalities typically associated with inspace operations are assembly (including maintenance of preexisting structures) and inspection. It is important to note that although treated separately, the assembly and inspection functionalities share common subtasks, such as navigation and maneuvering. Planetary surface exploration (PSE) functionalities include surface mobility (which includes navigation and collision avoidance), instrument deployment and sample manipulation, science planning and execution, and human exploration assistance. Even though space robotics applications are categorized into these two mission classes, generalization of the functionalities listed above allows for other mission profiles such as terrestrial assembly or inspace science. The research presented in the subsequent chapters is primarily applicable for inspace operations; therefore those will be the primary focus here. The area of inspace robotics operation is far less mature than robotic planetary surface exploration systems, both in terms of real world implementations and theoretical research. Real world implementations of planetary surface exploration have been significantly more numerous, and more successful, than inspace robotics implementations with examples of teleoperated surface exploration robots dating back to the 1970's with the Soviet Lunakhod Rover.3 The sophistication of more recent PSE system successes is evidenced by the Mars Exploration Rovers4 developed at NASA's Jet Propulsion Laboratory (JPL). In 1997 the Sojourner, although not designed for advanced capabilities, was the first Mars rover to successfully demonstrate the ability to teleoperate a vehicle on Mars. Most recently the two more sophisticated Mars rovers, Spirit and Opportunity, have demonstrated effective teleoperational capabilities in surface mobility and terrestrial science. The raison d'etre that led to this disparity between PSE and ISO implementations is twofold. First, the theory and experience of terrestrial robotics are directly applicable to PSE scenarios. Volumes of research since the 1950's have been devoted to the kinematics and control of terrestrial robot manipulators with achievements in autonomous navigation, perception, and mobility. Therefore, implementation of PSE robots spring boarded off the terrestrial robotics knowledge base and needed only to address the technical issues associated with high degrees of mobility, science planning, perception, and exploration capabilities, the harshness of space and space transport, and the high latency communication links implicit in space teleoperation. Conversely, with the exception of autonomous rendezvous and docking (ARD) investigations5 that date back to October 1967 in Russia6 and preliminary investigation of inspace proximity operations,7'8 it has only been in the last decade that autonomous, freeflying robotics applications have appeared in the literature.9'10 Therefore, both the underlying theory and its implementation need much further investigation in order to achieve the maturity of PSE robotics systems to date. Second, most terrestrial and PSE agents are developed without consideration of the system dynamics: they are typically kinematically driven, and the cases that do consider some dynamics do so only in the sense of constraints on the motion, e.g., nonholonomic constraints. For ISO applications, consideration of the system dynamics is essential. This fundamental difference is illustrated by the local operator in each case. In terrestrial or PSE applications, the local operator may be zero. This means that under zero command, the robot remains stationary until it is commanded to move. The local operator in ISO systems, however, is in general not zero. Instead, the local operator is the system dynamics, meaning under zero command, an ISO system will not remain stationary but move according to the topology shaped by the system dynamics. Therefore, consideration of the dynamics of ISO systems is imperative. The sections that follow outline the state of the art for each of the functionalities of ISO systems for both current implementations and experimental test beds. The expectations for the sophistication of these functionalities over the next decade are also discussed. State of the Art for InSpace Operations (ISO) Functionalities In this section both current state of the art and experimental test bed sophistication are assessed for assembly and inspection functionalities of ISO systems. The relatively few implementations of ISO systems to date are emphasized while noting the progress possible from the experimental investigations. The current state of assembly functionality, with respect to real world implementations, is restricted to the Shuttle Remote Manipulator System (SRMS) and the Space Station Remote Manipulator System (SSRMS) aboard the International Space Station (ISS). Neither of the current systems possesses any autonomy as they are solely teleoperated onboard or from the ground. However, to date there have been onorbit demonstrations of assembly subsystems such as the Robot Technology Experiment1 (ROTEX), a German experiment flown by NASA, and Engineering Test Satellite9 (ETS VII), flown by the Japanese Aerospace Exploration Agency (JAXA) (formerly National Space Development Agency of Japan (NASDA)). The ROTEX was a robotic arm (see Figure 1) that flew in 1993 on Columbia as part of STS55, and successfully completed multiple tasks that include replacement of a simulated Orbital Replacement Unit (ORU) and capture of a freeflying object via onboard and ground teleoperation and autonomous scripts. By accomplishing tasks from autonomous scripts during the experiment, ROTEX became the first autonomous space robotic system. The ETSVII mission was composed of a passive target satellite, Orihime, and a chaser satellite, Hikoboshi, with a robotic manipulator arm (see Figure 2). ETSVII successfully demonstrated cooperative control of a robotic arm and satellite attitude, and simple examples of visual inspection, equipment exchange, refueling, and handling of a satellite. Figure 1. Model of the ROTEX Figure 2. Artist's rendition of ETSVII manipulator components The state of the art in assembly robotics, implemented via experimental test beds only, is the Skyworker12 robot at Carnegie Mellon University (see Figure 3), NASA's Robonaut13 (see Figure 4), and Ranger14 at the University of Maryland (see Figure 5). Skyworker is an 11 DOF robot that walks across the structure it is assembling to mate new components to the existing structure. The current prototype allows for highlevel command inputs that are then parsed and implemented onboard as motion control commands. Robonaut is a collaborative effort between DARPA and NASA aimed at developing a humanoid robot capable of meeting the increasing requirements for extravehicular activity (EVA). Robonaut, composed of two dexterous arms and two five fingered hands with teleoperational and autonomous capabilities, has already demonstrated assembly of complicated EVA electrical connectors and delicate capabilities such as soldering. Ranger is a teleoperated robot at the University of Maryland that completes assembly, maintenance, and human EVA assistance tasks in a neutral buoyancy tank. Ranger has demonstrated robotic replacement of an Orbital Replacement Unit (ORU), complete endtoend electrical connector mate/demate, and twoarm coordinated control. Figure 3. CMU's Skyworker prototype Figure 4. NASA's Robonaut Figure 5. University of Maryland's Ranger According to "Space Robotics Technology Assessment Report""1 published by NASA in December of 2002, the expectation for inspace assembly under nominal research efforts is: "Robots that can autonomously mate components and do fine assembly, including making connections under careful human supervision."16 It is perceived, given the current state of the art of robots like Robonaut, that robots will possess the mechanical capabilities equivalent to a spacesuited human, but barring a breakthrough in communications architectures, more aggressive expectation of space robotic assembly cannot be met for teleoperation, due to current low bandwidth/high latency communication. On the other hand, automating Robonaut and other highly dexterous robots is possible, but to make highly dexterous robots effective under autonomous operation, a system level design needs to be considered that designs the small components specifically for assembly by autonomous robotic systems. This requires significant redesign in current infrastructure, which was perceived to be financially impractical at the time of the "Space Robotics Technology Assessment Report." However, the new NASA mandate outlines strategic allocation of funding to achieve its aggressive goal set, and as a result it now may be possible to alter the current infrastructure with the intent to simplify automated robotic assembly. Inspace inspection consists of examining space structures for anomaly detection. The robots performing the inspection can be either freeflyers or manipulators capable of performing the subtasks of moving, such as navigation and obstacle avoidance, to examine the entire exterior of the structure and anomaly detection via sensor interpretation. Currently, there are no operational space inspection systems. However, there have been onorbit investigations of subsystems, such as AERCam Sprinto1 (see Fig. 6), which flew aboard Columbia (STS87) in 1997. AERCam is a teleoperated EVA camera whose purpose was to merely investigate the feasibility of an autonomous, free flying EVA camera system. The state of the art in the area of inspection is the mini AERCam17 (see Fig. 7), Orbiter Boom Sensor System, 18 OBSS (see Fig. 8), and the Supplemental Camera and Maneuvering Platform, 19 SCAMP (see Fig. 9), at the University of Maryland. The mini AERCam project is a second generation version of the Sprint aimed at adding more complex capabilities while reducing the overall size of the prototype. A nanosatellite class spherical freeflyer, the mini AERCam is only 7.5 inches in diameter and weighs a mere 10 pounds. Even though mini AERCam cutting edge hardware is not flight certified, making it far from implementable, the results obtained from orbital simulations and 5 DOF experimental test bed validations are invaluable to advancing the field of autonomous inspection agents. The Orbiter Boom Sensor System (OBSS) is a manipulator based concept for inspection of the thermal protection surface (TPS) of the Shuttle. Consisting of a camera and a laser range sensor mounted on a 50' boom, the OBSS attaches directly to the SRMS, providing teleoperational inspection of nearly 75% of the Shuttle's TPS. In the Space Systems Laboratory at the University of Maryland, the Supplemental Camera and Maneuvering Platform (SCAMP) allows for investigations in freeflying camera applications. Operated in a neutral buoyancy tank, SCAMP provides near microgravity conditions for research into freeflying applications. SCAMP has demonstrated effective stereo video data interface and 3D navigation in the neutral buoyancy test bed. Figure 6. AERCam Sprint Figure 7. Mini AERCam Figure 8. Artist's depiction of the Orbiter Boom Sensor System ; " Figure 9. The SCAMP in the neutral buoyancy tank For inspace inspection, NASA perceives the expectation for the next decade to be summarized as "autonomous robotic inspection of some of the exterior surfaces with sensory data filtered for potential anomaly before being stored or sent."20 With more aggressive research effort, autonomous inspection and anomaly detection of most exterior surfaces of a target are realizable. The addition of autonomous capabilities to ISO systems is a timely issue. The new NASA mandate2 sets aggressive goals for space robotics missions, and autonomous capabilities for robotic systems are vital to the achievement of these goals. Furthermore, the international community, initiated by Deutsches Zentrum fur Luft und Raumfahrt (DLR), also known as the German Aerospace Center, has been proactive in forming a community whose sole focus is onorbit servicing (OOS). In their own words, the spirit of the OOS initiative "aims at developing a sound understanding of OOS and its issues, and at creating a dedicated community on global level"21 with the overall objective being "to develop a structured longterm approach of how to deal with OOS by providing high quality information and latest insights to this new, complex and innovative field of space activities, and to initiate new collaborations."21 OOS has organized two workshops, in 2002 and 2004, dedicated solely to all aspects of onorbit servicing, and numerous sessions for onorbit servicing in other conferences worldwide, e.g., the 8th International Symposium on Artificial Intelligence, Robotics and Automation in Space (Isairas 2005), providing further evidence to the point that ISO systems development is a timely, global issue. With the state of the art in ISO systems, as outlined above, reflecting few autonomous capabilities, and the absence of a mature knowledge base in freeflying robotics theory, adding autonomy to ISO systems is a monumental task as evidenced by the requests for information22,23 (RFI) issued from NASA in the months directly following the issuing of the new mandate, the recent formation of the Institute for Space Robotics, and the OOS initiative. To this end, this dissertation explores strategies for on board guidance and develops a guidance algorithm for proximity operations for autonomous ISO robotics systems in an attempt to improve the current state of the art of autonomous, inspace robotic systems. Research Scope For most ISO missions (assembly, repair, refueling) the target location and its orbit parameters are known a priori, with the exception being maneuvering noncooperative targets which we will not consider in the present work. For instance, in the case of an assembly mission, the location of the construction site and the parts depot are known parameters from mission planning. Utilizing this known information, we propose a novel approach to collision free operations of the spacecraft based on artificial potential function (APF) guidance. In order to utilize APF based guidance schemes, two prominent deficiencies inherent to this method need conquering: premature termination due to local minima in the operational space, and suboptimal performance. It is widely known that utilization of harmonic potential functions eliminates the introduction of local minima in the operational space; thus we employ harmonic functions to eliminate the introduction of local minima. In order to overcome the suboptimal performance, a sum of harmonic function primitives (simple functions such as sinks and vortices in the hydrodynamic analog) is employed to shape the velocity field of the operational space to harness the underlying dynamics of the system (i.e., orbit mechanics). Furthermore, given the location and orbit parameters of the target site, a random heuristic search algorithm can be employed to optimize the velocity field shaping for maximizing the performance gains achieved from harnessing these underlying dynamics. This dissertation presents this novel approach for proximity operations of multiple space assets operating in relatively unstructured environments. Dissertation Outline This dissertation is organized as follows. Chapter 1 provides an introduction to the current trends in the space industry, details the state of the art for ISO systems in both real world implementations and test bed environments, and concludes with the research scope that precisely defines the problem and solution presented in this work. Chapter 2 explores the orbit mechanics based approaches to trajectory generation for spacecraft, and illustrates the deficiencies that eliminate these algorithms for the specified application considered in this dissertation. In Chapter 3 the motion planning approaches from the robotics field are surveyed as alternatives to the classical orbit mechanics based approaches. From the robotics literature, the artificial potential function (APF) guidance approach is singled out as a promising solution approach for the problem considered here, and examples of APF guidance algorithms applied to the classical orbital rendezvous with and without obstacles are presented. Chapter 4 develops the novel alterations that overcome the two predominant deficiencies of APF algorithms that would prohibit their application for ISO missions. The first part of this chapter develops an original, truly general dynamic obstacle avoidance algorithm that is a critical subsystem for the trajectory generator for ISO missions. The second part of Chapter 4 deals with the mitigation of local minima and suboptimal performance by utilization of harmonic function primitives. The chapter is 13 concluded with the application of the APF guidance scheme to a representative ISO mission scenario. Eventually the spacecraft will not be considered to be a point mass and attitude considerations must be included. In Chapter 5 the classical APF framework is utilized to develop an attitude guidance algorithm that can be shown to be a subclass of controllers that can be derived using a Lyapunov analysis. Chapter 6 will present conclusions and future work. CHAPTER 2 TRAJECTORY DETERMINATION: ORBIT MECHANICS APPROACH The general mission scenario for all autonomous inspace operations (ISO) systems requires translating from an initial position to a goal position where the system fulfills its designed functional purpose, e.g., assembling a structure. It is assumed that the initial position of the system is either a parking orbit provided by the launch bus or an orbiting space platform (e.g., ISS), both of which imply the need for midcourse and endgame trajectory design. Midcourse trajectory determination refers to the transfer from initial location to close proximity, and endgame trajectory determination refers to the proximity maneuvering. In the sections that follow, classical orbit mechanics theory for the midcourse and endgame trajectory determination phases will be presented along with more recently developed methods based on orbit mechanics. Numerical examples illustrating the key aspects are included, and the limitations on each for autonomous ISO applications are illustrated. Midcourse Trajectory Determination The most general framework for midcourse trajectory determination is the nbody problem of astrodynamics (see Refs. 2427 or any text on Astrodynamics or Celestial Mechanics). The solution of which requires 6n integrals, of which there are ten known,24 27 and a general closed form solution is not obtainable. However, the artificial satellites are well within the sphere of influence (Ref 24, Chapter 7) of the primary central body, and the mass of the satellites is negligible with respect to the mass of the central body. Therefore, the trajectory of any robot in this context can be approximated effectively by the solution of a twobody problem, assuming the position of the goal is associated with an orbiting target. P: Figure 10. Schematic of Lambert's Problem Restricting consideration to twobody problems, the determination of the orbital trajectory given two points in space and a time of flight is known as the twobody orbital twopoint boundary value problem24 (see Figure 10). The first strides in the solution of this boundary value problem came from Johann Heinrich Lambert24'25 (17281779). Lambert's Theorem states that the time (t, t) required to traverse an elliptical arc between specified endpoints, P1 and P2, depends only on the semimajor axis a and two geometric properties of the space triangle FPP2, the chord c and the sum of the radii r, + r Since Lambert arrived at this conclusion from strictly geometric means, the mathematical statement of Lambert's Theorem, Lambert's Equation, actually comes from JosephLouis Lagrange24 (17361813) and is given by Eq. (1) in terms of Lagrange's variables a and/f. In order to show that Eq. (1) indeed proves Lambert's Theorem, a and / can be expressed as Eq. (2) and (3), where s = (r, + + c) is the semi 2 perimeter of the space triangle. Hence, the transfer time is a function of the semimajor axis, the radii sum r, + r, and the chord c. Obtaining solutions of Eq. (1) became known as Lambert's problem. Historically its solution was essential in determining the orbits of celestial bodies from observation, and today it is applicable for artificial satellite intercept and rendezvous problems. 3 (t 2 ,t) =a2[a,8 (sin ca sin)] (1) sin (2) 2 2a (3) sin K8 CJ (3) It is important to note that the solution of Lambert's Problem only defines the geometry of the intercept problem between points P1 and P2. To effectively execute the trajectory specified by the solution geometry, the terminal velocity vectors need to be calculated in order to obtain the necessary guidance command to initiate the transfer. Therefore, determination of the maneuver command is the essential quantity. Lambert's Problem The determination of an orbit from a specified transfer time (t2 t) connecting two position vectors, ri and r2, is called Lambert's Problem2427 (see Figure 10). A variety of techniques for the solution of Lambert's Problem have been developed over the years24,28,29 with each being classified by the particular form of the transfertime equation and the independent variable chosen in an iterative solution technique. Classically two transfertime equations24 were used, one from Lagrange's equations used in his proof of Lambert's Theorem and the other derived from Gauss' work. Each of these methods has their advantages and disadvantages, and a third transfertime equation that utilizes these advantages while mitigating the disadvantages can be developed as in Battin.24 Battin's combined equations technique is utilized for the present work. For a complete history and development of Lambert's Problem and the twobody orbital boundary value problem see Battin.24 Before attempting to solve Lambert's Problem it is instructive to consider inherent properties of Eq. (1) that will influence the techniques employed for its solution. Lambert's Equation is: 1. Transcendental equation 2. Double valued function of a (allows for conjugate orbits with same semi major axis) s 3. Derivative with respect to a is infinite for the minimum energy orbit a =  2 (the minimum energy orbit corresponds to the smallest possible value of semimajor axis a) 4. Four possible solution for each case due to quadrant ambiguities in a and/7 Property 1 requires the use of an iterative solution approach to find the semimajor axisa. From properties 2 and 3 it is obvious that the semimajor axis is not amenable to iterative solution techniques. Furthermore, Property 4 mandates the use of a method for selecting the unique solution from the multiple solution possibilities. These difficulties can be mitigated by utilizing a transformation of Eq. (1) by Battin24 that replaces the semimajor axis with an independent variable that produces monotonically decreasing and single valued solution curves, which is more appropriate for an iterative solution approach. Furthermore, utilization of Prussing's geometric interpretation25 provides the framework necessary to determine a unique solution. Battin's combined equations technique24 will now be presented. This procedure first transforms the independent variable in Lagrange's Equation to one with properties amenable to iterative solution techniques. Then, utilizing another form of the transfer time equation attributed to Carl Friedrich Gauss (17771855), Battin improves the computational efficiency by reducing the number of hypergeometric function evaluations from two to one, however the independent variable in this transfertime equation lacks the attractive properties obtained previously. By developing a relationship between the two independent variables, Battin combines the two transfertime equations in order to obtain the combined equation form utilized in the current work. Battin's combined equations transfertime equation is given by Eqs. (4a) through (4f), where F [a, f; y; x] is the hypergeometric function. This form of the transfertime equation is devoid of the difficulties encountered from using the semimajor axis as the independent variable, which renders it readily adaptable to an iterative solution method, such as Newton's Method. The parameter A is a function of the known geometry and y is an intermediate variable expressible in terms ofx as Eq. (4c). Furthermore, ris essentially a function of x only since y can be given by Eq. (4c), and S, is a function ofx only, which implies that Q = Q(x) only. Therefore Eq. (4a) is a function of the independent variable x only. It is also important to point out that the solution obtained for x is unique, a property not apparent by inspection. The reason for the uniqueness can be attributed to the use of variable x andy, which can be given by alternate expressions 1 1 x = cosa and y = cos/ Therefore, ambiguities arising from the trigonometric 2 2 functions involving a and / are avoided. I(t tl)= 3Q +4A)7 (4a) V am c= cos (4b) s 2 y 1= 2 x2) (4c) S= y Ax (4d) S = (1 A 77x) (4e) 2 Q =F 3,1; 5; S (4f) 3 2( Maneuver Command Determination In order to calculate the maneuver command the solution x from the Battin algorithm must be converted to geometric properties a, a, and f. The equations utilized S for this purpose are given by Eqs. (2), (3), and (5), where a is the value of the semi 2 major axis for the minimum energy transfer between Pi and P2. From the inverse trigonometric operations the conversion always yields the principal values ac and /0, characterized by 0 < /0 < o a < arising from quadrant ambiguities from inverse trigonometric operations. In order to obtain the correct solution from the four possibilities, Prussing geometric interpretation25 is employed. x2 am (5) a Pressing geometric interpretation involves transforming the geometry of Figure 10 to a rectilinear ellipse, thereby introducing geometric significance to the angles a and . The results of this procedure provides a framework for the logical determination of the appropriate values of a and / from the principle values oa and /0. The assignment procedure is given below by Eqs. (6) and (7), where (t2 t1) is the transfer time, tm is the transfer time on the minimum energy transfer, and 0 is the transfer angle. a {2 a ,( = a 6, o 2 ) \ 0,O Once a and / have been obtained from the above procedure the terminal velocity vector v, at r1 (andv2 at r2 for a rendezvous) is calculated. To do so, introduce a skewed set of unit vectors, defined in Eqs. (8a) through (8c), where i~, i = 1, 2, are unit vectors in the direction of r,, and ui is a unit vector parallel to the chord. It can be shown25 that in this coordinate bases v, has the elegantly simple expression given by Eq. (9), with A and B given by Eqs. (10) and (11). Therefore, the final objective of obtaining the maneuver command necessary to initiate the trajectory between the start and goal configuration can be calculated by taking the vector difference of the satellite's current velocity with v,. r 1 (sa) '1 u2 = 2 (8b) r2 S= 2((8c) c EI = (B+ A) fc ( A)fi A = cot (10) (4a) 2 B = C cot (11) 4a 2 Numerical Example: Midcourse To illustrate the midcourse trajectory determination and maneuver command algorithms, a numerical example is presented. The problem is as follows: given position vectors r = [9567 0 0] and r2 =[37394 10783 0] (units are kilometers) coordinatized in the Perifocal (PQW) frame (xaxis pointed in the direction of the eccentricity vector) and a time of flight equal to four hours, determine the solution trajectory and the maneuver command. The initial position was chosen to be the periapse which allows for simple analytical verification of the terminal velocity vector produced by the midcourse maneuver command algorithm. Using algorithms based on the theory of the previous sections and developed by the author the solution is obtained. The iterative solution for x in Lambert's Problem using Newton's method required 19 iterations to converge with a relative error less than 108. The terminal velocity at P] was found as v, = 8.73 xlO10 8.228 01T kilometers per second. The CPU runtime on a dual 2.4 GHz Xeon processor workstation for the midcourse maneuver command determination algorithm was approximately 0.15 0.05 seconds. The terminal velocity v, produces by the algorithm can be verified by checking i (I+ e) against the analytical expression for the velocity at periapse passage v = (1e Performing this calculation produces the analytical terminal velocity vector as Performing this calculation produces the analytical terminal velocity vector as p = [0 8.228 0] km/s, which shows good agreement with the results from the algorithm. Figure 11 depicts the transfer trajectory, obtained from a Keplerian motion simulator using the initial conditions given by the position and terminal velocity vector solutions. It is important to note a caveat for the numerical example above. The position vectors were given in the Perifocal (PQW) frame, which is an orbit based coordinate frame. However, it is unlikely that the position information will be known in the PQW frame, instead position information is more likely to be obtained in the EarthCenter Inertial (ECI) frame or the LocalVertical LocalHorizontal (LVLH) frame. The implication of this is, unlike the example above, the initial value of the true anomaly is unknown and cannot be determined by inspection from the position vectors. For instance, if the position vectors in the example above were given in the ECI frame (with Q = 10,i = 28.5, co = 45) the position vector would be r1 = [5629.4 7694.7 792.9]' and ri = [14080.1 35350.3 8166.2] km, which does not exhibit any evidence that the true anomaly is zero at P]. Therefore, in most cases only the difference in true anomaly is known from geometry and the elegant calculations used to verify the terminal velocity vector here are not applicable. Endgame Trajectory Determination Endgame conditions for this work will be defined such that the dynamics of relative motion between two bodies in orbit can be accurately approximated by linear orbit theory. The limit on the applicability of linear orbit theory is dependent on both the desired level of fidelity and the class of orbit operation, e.g., LowEarth orbit (LEO) or Geostationary orbit (GEO). For an indepth analysis on the applicability of linear orbit theory see Refs. 30 and 31. For the analysis that follows it is assumed that endgame conditions are met. x 104 2 Solution STr ajectory : p r,= =9567 Anm 1: 1.5 0. 1 ........ . : =392 8160 o x [km] Figure 11. Transfer trajectory from midcourse trajectory determination algorithm Chaser Figure 12. Depiction of orbital rendezvous The basic scenario for the development of linear orbit theory is depicted in Figure 12, where both the Target and Chaser are assumed to be in twobody motion about the x 10 central body, denoted by M. The most familiar equations from linear orbit theory and the ones employed in the present work are the ClohessyWiltshire (CW) equations,32 originally developed for the first orbital rendezvous flight demonstration programs. The CW equations are given in state transition matrix form by Eqs. (12) (14), where x, y, z (z out of the page) are the coordinates of the Chaser with respect to the Target in the local frame depicted in Figure 12 and n is the orbital mean motion of the Target. It is vitally important to understand the implicit assumption of the CW equations that ultimately limits their applicability and fidelity. The reference (Target) orbit is assumed to be circular. Therefore, for any scenario where the reference orbit is not circular errors in addition to the ones from linearization are introduced, thus, further degrading the fidelity of the model. r() = ()r(0) (12) r (t) =x(t) y(t) z(t) x(t) y(t) z (t) (13) s 2 43cosnt 0 0 (cosnt) 0 n n 2 4 sin nt 3nt 6(sinnt nt) 1 0 (1cosnt) 0 n n S(t) 0 0 cost 0 0 (14) n 3n sinnt 0 0 cost 2 sinnt 0 6n(1cosnt) 0 0 2 sinnt 4cosnt3 0 0 0 n sinnt 0 0 cost CW Maneuver Command Determination The primary advantage of the CW equations is the computational simplicity of determining the transfer velocity necessary for performing a maneuver. The objective of the maneuver is to bring the Chaser with initial conditions 8r(0) and 8r(0) to the origin of the CW frame (depicted in Figure 12) given a time of flight t*. Partitioning the state transition matrix D into four 3 x 3 matrices as given by Eq. (15), the top partition leads to the equation for the relative position r(t) = M(t)dr(0)+ N(t)d(0) Substitute t =t*and dr(t*)= 0, giving the transfer velocity r,0 the expression shown by Eq. (16). M(t) (t) N(t) (15) (t ) =( (15) ( S(t) T(t)J + =N(t*) M(t*)dr(O) (16) Numerical Example: Endgame Continuing the maneuver begun in the numerical example of the midcourse section, the problem scenario is as follows. The robot completes the midcourse trajectory to find the goal in a circular orbit with radius rG = 38916.7 km. The discrepancy in the positions of the goal and robot at the end of the midcourse maneuver can be viewed as errors incurred from unmodeled perturbations such as drag, oblateness, or solar radiation pressure effects. The objective of this example is to use the endgame algorithm outlined in the previous section to generate the maneuver command that will bring the Chaser to the Target's position with a time of flight arbitrarily specified as onequarter of period of the Target orbit. All vector quantities are coordinatized in an Rbar Vbar frame attached to the Target, with x in the direction of Rbar and y in the direction of Vbar. The endgame maneuver command generated by the author's algorithm is given by Av= [0.0001 1.4028 0] km/s with a CPU time less than a few milliseconds on the workstation used previously. The maneuver was simulated, as depicted by Figure 13, and the absolute value of the errors in position at the end of the maneuver were obtained as I = 1.059x 10 9 1.206 x 109 0 which illustrates the convergence of the robot to the goal position in the specified time. 0.05 Ini ial Endgame Goal Position 0.05 SEndgame 0.1 trajectory 0,2 0 0. 02 0.3 0.4 05 0.6 0.7 0.8 0.9 1 x [kml Figure 13. Endgame maneuver trajectory in the xyplane (Rbar Vbar plane) Alternative Orbit Mechanics Based Methods Having illustrated the classical methods of orbit mechanics for trajectory determination, we turn to a more recent algorithm33 based on an adaptation of glideslope guidance that has been utilized for rendezvous and proximity operations for the space shuttle and other vehicles with astronauts in the loop. To date this is the most widely used approach from a real world implementation perspective. The glideslope approach utilizes the linearized equations of motion Eqs. (12) (14) but completes the maneuver with Nthruster firings. The algorithm divides the rangeto go into N equally spaced segments along the chord connecting the initial and final relative positions. At equally spaced time intervals an impulse is calculated using the CW equations that drive the current position to, the next intermediate position rm. An example of the algorithm is shown in Figure 14. x [km] Figure 14. Glideslope solution for 5 impulse rendezvous. The computational simplicity of the glideslope algorithm lends itself readily to application to ISO mission. However, this computational simplicity is sacrificed for situations where obstacle avoidance is necessary. In order to use a glideslope approach with obstacle avoidance, a chord must be drawn between the centers of mass of all the spacecraft in the operating environment. Therefore, (n1)! chords, each with Nimpulse to be calculated are necessary. For instant for a five spacecraft mission, requires 24N online impulse calculations! For this reason alone the glideslope approach is an unacceptable alternative for the ISO missions considered in this dissertation. Conclusions The existence of two trajectory determination algorithms, one for midcourse and one for endgame, implies that a multitier architecture is necessary. Therefore, not only are two algorithms necessary but the added complexity of switching algorithms to determine when the linear control is valid is introduced. To mitigate a mutitier control structure, the midcourse maneuver commands could be utilized for both midcourse and endgame scenarios. However, the midcourse maneuver command determination is a computationally intensive algorithm that requires an iterative solution, which does not lend itself to onorbit implementation due to the limited computational power that exists for current space systems. Another alternative is to modify the linear orbit theory in order to expand its range of applicability to midcourse trajectory determination as well. This subject has received intense research,3437 including investigations by the author, with minimal achievements in expanding the range of applicability. An additional disadvantage of classical orbit mechanics theory was obvious to Stern and Fowler30 twenty years ago; classical orbit mechanics theory does not admit the possibility of the incorporation of physical trajectory constraints, introduced by either goal structure or obstacles in the path. Since the ability to dynamically incorporate obstacles and trajectory constraints is vital for all autonomous ISO systems, and given the additional disadvantages mentioned above, other sources of trajectory planning theory for autonomous ISO systems must be explored. Chapter 3 will survey the robotics literature for motion planning algorithms and will motivate and develop artificial potential function guidance algorithms for autonomous ISO applications. CHAPTER 3 ARTIFICIAL POTENTIAL FUNCTION POSITION GUIDANCE In the preceding chapter it was demonstrated that trajectory determination for autonomous inspace operation (ISO) systems cannot be easily handled by classical orbit mechanics theory. A logical avenue to explore for trajectory determination is the robotics literature on robot motion planning (MP). A brief outline of the robot MP problem is discussed in the first section of this chapter with an emphasis on motion planning approaches. Given that spaceborn robotics applications require motion planning subject to properties distinct from classical robot MP problems, for instance the difference in local operator mentioned in Chapter 1, only some of the common motion planning approaches are feasible for inspace applications. Based on the unique needs of inspace applications, a potential field approach to MP is perceived to be most appropriate for autonomous ISO systems, and this chapter concludes with a presentation of classical artificial potential field (APF) guidance algorithm for ISO applications. Robot Motion Planning In order to precisely define the robot MP problem, it is first necessary to define two terms, configuration and configuration space (Cspace). A configuration of a robot is defined as the set of independent parameters, degrees of freedom (DOF), necessary to specify the location of every point of that robot. Cspace, a concept introduced in an influential paper by LozanoPerez and Welsey38 in 1979, is defined as the set of all possible configurations of a robot, i.e., the Cspace represents all possible motions of a robot. Essentially all MP problems are equivalent when cast in the Cspace, and can be stated as follows. Find a connected sequence of point in the Cspace between an initial configuration and goal configuration. Configurations that result in collision with the environment, known as configuration obstacles,39 must be computed, and the connected sequence of points must not contain any points associated with configuration obstacles. Computation of the configuration obstacles is difficult39'40 and the dimension of the Cspace is equivalent to the DOF of the system, both of which makes motion planning challenging for robots with multiple degrees of freedom. Classically the MP problem has been broken down into the steps below. Concentration on the approaches in the literature for Steps 3 and 4 are the focus here for determining an applicable approach for ISO systems. For a complete investigation on methodologies for the steps in the MP problem see the excellent surveys in Refs. 3943. 1. Parameterize the configuration of the robot(s) 2. Choose a representation scheme for robot(s) and the environment 3. Select a motion planning approach 4. Select a search method to find feasible path through environment 5. Optimize the solution path based on various constraints (e.g., smoothness) Citing the landmark survey by Hwang and Ahuja,39 the most common motion planning approaches found in the literature are; cell decomposition, mathematical programming, skeletons, and potential field methods. Cell decomposition techniques involve computation of the entire Cspace and all configuration obstacles in addition to requiring large memory allocations. Therefore, given that onorbit computational power is limited, this approach is not feasible for ISO systems. Mathematical programming methods represent the requirements for obstacle avoidance as a set of inequalities on the configuration parameters. The motion planning problem is then formulated as an optimization problem that finds a solution path between start and goal configurations based on minimization of a scalar constraint. This motion planning approach is essentially a nonlinear optimization problem subject to multiple inequality constraints, which generally requires a numerical method for its solution. Again this process is computationally intensive and therefore does not lend itself to onorbit implementation. Skeletons, otherwise known as roadmap or highway approaches, contract the environment into a network of 1D lines and the motion planning problem becomes a graph search problem restricted to the network. Common examples of the skeleton method are Voronoi diagrams and subgoal networks. The skeleton methods have low memory requirements, and if an accurate representation of the environment and obstacles is known a priori, this method is feasible for ISO applications. In fact a Voronoi diagram approach is currently being investigated for the AERCam project.44 Potential field methods for motion planning are the simplest, computationally speaking, of all the motion planning approaches presented here. This method creates a scalar function called the potential that has the following properties. * Potential is at a minimum when the robot is at the goal configuration * Potential is high on the surface of obstacles * Everywhere in the environment away from the obstacles, the potential slopes toward the goal configuration Typically the potential field method is combined with a specified local operator, such as "gostraight" and the path taken to the goal is recorded. Potential field methods also allow for simple consideration of nonpoint robots by evaluating the potential at representative points on the robot. Additionally, no explicit Cspace computation is necessary. It is the opinion of the author that the utilization of an artificial potential field motion planning approach with the local operator being the dynamics of motion is the most effective motion planning approach for ISO applications. Artificial Potential Function Position Guidance Artificial potential function (APF) guidance for satellite systems4549 has appeared in the literature recently. However, the application of realtime proximity operations of an autonomous ISO system has yet to be investigated. In this section the robot will be considered to be a pointmass moving in three dimensions. This restriction will be removed for attitude considerations in Chapter 5. The objective of the APF position guidance algorithm is to drive the robot to the origin of the state space through a series of impulsive maneuvers. It is instructive at this point to state the equation of motion of the dynamical system considered in this section. Since this section is focused on the relative motion between a robot and target in a twobody orbit, the motion of each is governed by Eq. (17), where g (r) is a vector function describing the gravitational field and aD is the disturbing accelerations due to perturbations such as atmospheric drag or oblateness. It will be assumed that the perturbations are zero and the both satellites are in an inverse gravitational field given by g(r) = r3 for the present analysis. Therefore, the relative motion between the robot and target &r = rR r is given by the vector difference of the governing equation of the robot and the target, as given by Eq. (18). However, these are the equations of motion in a Newtonian reference frame. It is preferable to transform these equations into a local vertical local horizontal (LVLH) frame attached to the target. This is accomplished by the classical transformation given in Eq. (19) where primes denote differentiation with respect to the local coordinate frame and o) and c are the angular velocity and acceleration of the local frame, respectively. Thus, the equations of motion for the relative motion of the robot and target are given by Eq. (20). Furthermore, it is obvious that the origin of the LVLH frame 0r = 0 is indeed an equilibrium point of the dynamical system. r = g(r) + aD (17) iL = R i = g(R)g(rT) (18) r = r" + 2o x r'+ o' x r +o) x o) x r (19) r" = Pr + rT (2oxr'+co' xr+cox coxr) (20) rR rT The general theory for APF guidance is based on the 2nd method of Lyapunov5s which can be described as follows: Consider the dynamical system x = f (x, t) with x e 9 An equilibrium point of the system x, is such that f(xe,t) = 0. For simplicity it will be assumed that through a change of coordinates any equilibrium point can be shifted to the origin, and therefore, Lyapunov's theorem will be presented for equilibrium states at the origin only. Essentially, Lyapunov's 2nd Method is a means of assessing the stability of the equilibrium points of the dynamical system.50 The theorem states that if there is a scalar function V (x,t) such that Eqs. (21) (23) are satisfied, whereV V (,t is the derivative evaluated along the solution trajectories, then the origin of the state space is a globally attractive point. V(x,t) > Vx#0 (21) V(x_,t) > as IIX 4 (22) r t <0 VxO 0(23) The first two conditions are satisfied by choosing V to be an appropriate positive definite function of x (it will be assumed for now that the potential function is not an explicit function of time). The third condition in Lyapunov's 2nd theorem is enforced by an appropriate choice of control action.51 The total derivative of the potential function is given by V = (VV)T x, and the control is defined by Eq. (24) where 8 is such that the velocity after a control impulse x+ is directly opposite to the local gradient of the potential _x, = kVxV. This proper51 control constrains the motion of the robot after an impulse to be tangent to the local gradient of the potential. Defining the control in this manner ensures the derivative of the potential after a control impulse is negative definite, as given by Eq. (25), and therefore Eq. (23) is satisfied. Thus the origin is a globally attractive point under the control action specified. Control V< (24) \8, V>0 V'(x)= k(VV) (VV) <0 (25) A pseudocode for the APF guidance algorithm is given below, and an example of APF guidance for an unconstrained (no obstacles) rendezvous is given in the next section to demonstrate the implementation of the algorithm. Algorithm 1. Obtain the current relative position and velocity between robot and target 2. Check the switching condition V = (Vv)' x 3. If V > 0, give an impulsive velocity input such that x,+ = kVxV 4. If V < 0, coast under the influence of the dynamics of motion 5. Sample position and velocities at desired frequencies and apply control based on switching conditions 6. Repeat until terminal error bound is satisfied Example: Spacecraft Unconstrained Rendezvous Define the potential function V(x) = x Px, where P is assumed to be a positive 2  definite, diagonal matrix for simplicity and x = &r = [x y z]T are the coordinates of the relative position in the LVLH frame. This gives VxV = Px and the total derivative S= (px) x For the rendezvous, the relative velocities of the spacecraft x must also converge to zero. The convergence is actually guaranteed by the choice of control presented above. Consider the norm of the velocity after an impulse given by Eq. (26), thus _d > 0 as Ix > 0 and the rendezvous is guaranteed. _l = k[Px]2 (26) The numerical example considers the rendezvous between a robot in a circular orbit with radius rR = 6705 km and a target satellite also in a circular orbit with radius r, = 6700. The robot leads the target radius vector by a phase of 0.001 degrees. The APF guidance parameters used were k = 1 and P = I, where I equal to the 3 x 3 identity matrix. Figure 15 depicts the trajectory taken by the robot in the xyplane, where x and y are the coordinates of the LVLH frame in the orbit plane with x in the radial direction. The sharp corners on the trajectory mark where an impulsive control was applied. Figure 16 shows the impulses in the x and y directions, respectively, applied by the APF guidance algorithm as a function of time. In all a total Av = 27.10 m/s was necessary to bring the robot to within a meter of the target. The duration of the maneuver was 86.2 minutes. It is important to note that the APF guidance algorithm ensures asymptotic convergence to the origin, thus a majority of the duration of the maneuver came from closing the final meters. 41 4  1 0 1 2 3 4 Figure 15. Planar rendezvous in the CW frame under APF guidance A010 It 5011 IWoC Iwlo] OOfWj (1 Sa.a a w5 0 w0 I Figu 15 P t in the frm un A gu 0 500 O 1500 M2O 25M 300 3 4000 45 SOw .M f i S; Figure 16. Control efforts in x and y directions for the rendezvous Obstructions In proximity operating conditions, instances exist where certain regions of the state space need to be avoided; for example, solid obstructions in the flight path. Therefore, the potential function needs to be altered in order to ensure the robot avoids these forbidden areas. For the current work, all obstructions will be approximated by spheres (or circles in the planar case). This is a reasonable assumption for two reasons. Firstly, for now the objective of the guidance algorithm is to provide safe trajectories between a start to goal location, not to generate trajectories for close proximity flybys, yet. Therefore, an approximate map of objects is acceptable. Secondly, the big picture for autonomous ISO systems is perceived to start with a rough estimate of the working environment, with perhaps this estimated environment approximating objects initially as spheres and updating the map through sensor data obtained from flybys. The obstructions will be represented as regions of high artificial potential in the potential field in the form of Gaussian function as in Eq. (27). The parameters Vrk and ko define the width and skewness of the Gaussian function whose center is located atxk, and the matrix M defines the shape (since all obstacles are assumed spherical M = I). Representing obstacles by Gaussian functions ensures that controls remain bounded, no singularities are introduced into the potential field, and there are no local minima. It is instructive to note here the units of the width and skewness are [L2] which will be important due to the fact that all simulations are run in nondimensional units for completeness. 'k (xk)= vk exp {r (xk )M(x k)} (27) Once the APF for the obstacles has been selected, the parameters of the obstacle APF must be determined such that the robot cannot cross into the forbidden regions. At first glance it would appear that constraining the potential everywhere on the surface of the obstruction to be equal to the potential at the initial conditions, stated mathematically in Eq. (28), would be sufficient. However, this requires the width and skewness to be function ofx, which would degrade the most attractive property of the APF approach, simplicity. An alternative constraint such that the V(x,,,) 2 V(_x) still guarantees the robot cannot cross the obstruction surface, and constant values of y'k and o can be found to satisfy this constraint. V (x) V(X ,) V(xo) (28) Utilizing the assumption of spherical obstructions only, a simple analytical method for the determination of the obstruction APF parameters rk and ok can be obtained. The potential on the surface of the obstruction is given by Eq. (29). For spherical obstruction (s_,ur )xk surf k) = R where Rk is the radius of the spherical obstruction. Effectively there are two parameters to be chosen based on one constraint equation, thus only one independent parameter. The prescription for obtaining yk is to choose a value for o and vary x,, over all the points on the surface of the obstruction calculating a Vrk at each value via Eq. (30). Then looking at the resulting potential fields for the various ryk, select the ones that do not violate V(x_ ) > V(x0o). A simulation was implemented to investigate which values of x,,s produced Vyk such that the constraint holds, the results of which are illustrated by Figure 17 through Figure 20. V(x,,,) = x Px,,,f +k exp {crk1 ( x ) (x x) (29) Vk = exp R V(x) x Px ,,) (30) Figure 17 shows percent error between the potential at the initial conditions and the potential on the surface of the obstacle. The black line represents the percent error incurred when using the value of Vfk calculated using xs,, for which ixs is a maximum over all points on the surface, and the red line corresponds to percent error using the value of lVk computed from x,, such that _xJ is a minimum over all points on the surface. It is clear from this diagram that the latter value of xsuf must be used to compute Vrk in order for V(x,, ) > V(x,) to be satisfied. For spherical obstructions finding xs, whose norm is a minimum over all points on the surface x* is a simple exercise in geometry, and can be computed by Eq. (31) were x is the vector connecting the center of the obstruction to the origin. xsurf = 1I k (31) Figure 18 and Figure 19 depict contour plots of the potential in the xyplane using Vyk calculated by Eq. (31). Figure 18 shows the level of contour at the initial conditions as labeled. Figure 19 is a zoomed view of the previous figure in the region of the obstacle (for this case the obstacle center was given by k= [2 0 0] ). The blue doted line corresponds to the positions on the surface of the obstacle, and it is evident that the potential is greater than the potential at the initial conditions at all the points on the boundary. Figure 20 shows the percentage error between the potential at the initial conditions and potential on the surface of the obstruction foryfk as calculated by Eq. (31). Figure 17. Difference between IC potential and potential on surface of obstructions catws Mop of Pdn"W Fhotwf Obapmbmn W [MO.M X lnri Figure 18. Contour plot of the potential function Eirar Bthu V, Wd Rad PM Obtinrtn Pote X khl Figure 19. Zoomed view of contour plot Obstrueltan &uioe PoUtit Earr vs. Suface Poaion 12, I 2: I T Y Inrl 1 1 2 Figure 20. Stem plot of obstruction surface potential error 23 24 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 x [knl Figure 21. Trajectory traversed by robot for the pathconstrained rendezvous example Numerical Example: PathConstrained Rendezvous The initial conditions for the previous example will be used again in this example. The APF guidance algorithm parameters used were k = 0.25 and P = I3, the different value of k was selected to enhance the illustration of the behavior of the robot around the obstacle. The obstacle is a circle in the plane with center x = [2 0 0] and radius R = 0.5 km. The parameter rk = 0.01 which leads to ak k = 1.0855 x 105, both given in nondimensional units (the nondimensionalization employed the radius of the target as the length scale factor and the mean motion of the target as the time scale factor). Figure 21 shows the trajectory traversed by the robot during the constrained rendezvous maneuver. The duration of the maneuver was approximately 245 minutes and a total Av = 61.52 m/s to bring the robot to within a meter of the target. Figure 22 (a) and (b) and Figure 23 (a) 43 and (b) depict the impulses in the x and y directions, respectively, applied by the APF guidance algorithm as a function of time. The second plot in both figures is zoomed around the times the robot encounters the obstacle. x Io(a) "q3 Til[7 1400 1600 1800 2000 2200 2400 2600 Tiw I[s] Figure 22. Control action in the xdirection for the constrained rendezvous (a) entire simulation (b) control actions in vicinity of obstacle (0) (a) l1.01 ......... 0 50 100oo00 I'd '.1. 515 15oo 2000 15000 Time ls] Figure 23. Control action in the ydirection for the constrained rendezvous (a) entire simulation (b) control actions in vicinity of obstacle "` ~~~ Conclusions In this chapter the motivation for employing artificial potential function (APF) guidance for autonomous ISO applications was presented. Furthermore, position guidance based on APF was developed for unconstrained rendezvous and path constrained rendezvous, with a numerical example presented for each. Through the numerical examples it was illustrated that APF position guidance can be an effective tool for autonomous proximity operations, a vital element for all ISO systems. Even though the classical APF examples from above were shown to be effective, they suffer from two predominant drawbacks of APF methods: local minima causing premature termination of algorithm, and suboptimal performance. Furthermore, the absence of a truly general APFbased obstacle avoidance algorithm needs to be addressed. In the next chapter, the APF guidance framework will be augmented with novel alterations that both conquer these drawbacks and extrapolate the obstacle avoidance algorithm to the most general dynamic case. CHAPTER 4 ARTIFICIAL POTENTIAL FUNCTION GUIDANCE: A NEW APPROACH The earliest implementation of APF algorithms is attributed to Khatib and Mampey52 in 1978, who used the approach for force control applications. Since then, APF methods have been applied to a wide variety of research areas, most predominantly robot motion planning5355 and computer graphics.56 Recently, APF guidance has been utilized in onorbit applications such as autonomous rendezvous,7 position,5 and attitude59'60 guidance for inspection, maintenance, and assembly operations; and formation control. Despite the proliferation of APF methods in the robotic literature and their diverse utilization in robotics research, the two primary drawbacks of these methods, spurious local minima and suboptimal performance, have received sparse attention. Furthermore, in the vast body of literature on APF methods a void exists for algorithms capable of handling a general dynamic environment. The first section of this chapter will address the vacancy in the literature pertaining to APF based obstacle avoidance algorithms for truly general dynamic environments. Section two of this chapter will be charged with overcoming the two predominant shortcomings by utilizing a sum of harmonic function primitives to compose the attractive potential function. Dynamic Artificial Potential Function Obstacle Avoidance Algorithm Despite the numerous applications of APF methods, most of these instances deal with static environments (potentials). However, a few examples of published works5943 dealing with an APF approach in dynamic environments exist with the most prominent being the work of Ge and Cui.63 Nevertheless, these works introduce prohibitive restriction as to the nature of the motions of the target and obstacles (e.g., constant velocities for both). Furthermore, they only consider the situation were the robot is required to take obstacle avoidance measures for one obstacle at a time. To the author's knowledge, a truly general APF algorithm for dynamic environments without these restrictions has not been addressed. The dynamic artificial potential function guidance scheme employed in this section consists of two components: an attractive and a repulsive potential. The attractive potential is a familiar form seen in previous studies on artificial potential field applications. However, the repulsive potential has been augmented with a priority weighting scheme that implements simultaneous obstacle avoidance for all obstacles requiring such. This weighting scheme enables APF guidance to be extended to dynamic environments, with the caveats that the robot possesses sufficient actuation to complete the maneuver, and the refresh rate for updating environment data is sufficiently fast. Attractive Potential The attractive potential, Uatt, given by Eq. (1) includes both the position, r (t), and velocity, v(t) of the chaser and position and velocity of the target, rtar (t), vt (t) respectively. The constants m, n > 1 and a,, a are design parameters. U,, =a lr,, (t)r(t) m + a, IIv, (t)v(t) (1) Similarly to the classical APF approach, we define the virtual attractive force, Fatt, as the sum of the gradients of the attractive potential with respect to the position and velocity of the chaser, as given in Eq. (2). Evaluation of two gradient operations yields the expression for the virtual attractive force shown by Eq. (3), where nRT is a unit vector pointing from chaser (robot) to the target (goal) and n', is a unit vector in the direction of relative velocity between chaser and target. It is clear to see that the virtual attractive force goes to zero as the chaser simultaneously obtains the same position and velocity as the target. F tt = VpUa VU,, (2) Fa, = maCp rtr(t) r(t) lm l1 n + na, lcVt (t) v(t)l n, (3) Repulsive Potential In order to define the repulsive potential, parameters need to be introduced to quantify the necessity for obstacle avoidance as it pertains to each obstacle in the environment. Define the parameters given by Eqs. (4) (8) for each obstacle i = 1...n. With these parameters in mind, we define the repulsive potential as given by Eq. (9). Qualitatively, the potential is zero when the chaser is moving away from the obstacle, as indicated by 'vRo (t) < 0, while the potential is nonzero for situations when the chase is approaching an obstacle and is within that obstacles range of influence. The repulsive potential is not defined when minimum distance between chaser and the obstacle is less than pm, meaning a collision cannot be avoided. The total repulsive potential is obtained by summing each repulsive potential over all the obstacles, Urp = 'Ue . amax maximum acceleration of chaser (4) p, = min distance between robot and obstacle = r (t ' ro (bt) 'r (5) 'Pm distance under amax 'VRo 0 RO (6) 2amax 'PA influence range of obstacle= 'pO ('vR (t), T) (7) VRO(t)= v(t) 'vob (tT n RO, = ...n (8) 0, if 'p, > 'p0 or 'vo (t)< 0 'UI = { r 'Pm if 0 < 'p, 'pm < and 'vR (t) > 0 (9) SIP P. P, Po.'0) not defined, if 'vRO (t) > 0 and 'p, < 'Pm Just as in the attractive force case, the virtual repulsive force is defined as the sum of the gradients of the repulsive potential with respect to the position and velocity of the chaser. After performing the operations and simplifying, the virtual repulsive force, which invokes the necessary obstacle avoidance, is given by Eqs. (10) (12), where 'VRO is the magnitude of the component of the relative velocity of chaser with respect to the ith obstacle in the direction perpendicular to the line connecting the chaser and that obstacle. 0 ,if 'p, > 'po or 'vR (t) < 0 F rep = Frep, +'Frep2 if 0 <'p 'p < oand 'vRO (t) >0 (10) undefined if'p, <'pmand 'vRO (t) > 0 rep r, ep VRO(t) O (11) Samax 'Fep2 q p vRO (t) vRO "O (12) amax Ps Given the repulsive force definition above, the priority weighted scheme is elucidated. Introduce the priority index for the ith obstacle ('PI), calculated as shown by Eq. (13). Again, qualitatively the priority index is a measure of how close are the minimum distance under action of maximum acceleration and the shortest distance between the robot and obstacle. As any one 'PI approaches some predefined factor of safety, FOS, (one would prefer not to tempt fate by waiting until the PI 1 and the robot was forced to do a critical avoidance maneuver) the obstacle avoidance algorithm switches from the simultaneous mode to a single obstacle avoidance mode for the obstacle with critical valued PI. Note that we did not consider the case of simultaneous obstacles with critical valued PIs. This logic guarantees a collision free trajectory through the operational space. It is important to note at this juncture that the priority index implicitly depends on relative velocity along the chord connecting the two centers of mass of the robot and obstacle. This implicit dependence is crucial to the situation with nonconstant velocity targets and obstacles. S if 'ps p > po or 'vR < 0 'PI = p (13) Pm' otherwise (13) Under normal operating conditions, i.e., no 'PI> FOS, the total repulsive force felt by the robot is calculated as the PIweighted sum of the individual repulsive forces for each obstacle in the environment, as shown in Eq. (14). Frep = PI* 'Frep (14) Total Force Interactions The total virtual force acting on the chaser is the sum of the attractive and repulsive forces, Fo,,, = Fa + Fep Through Newton's 2nd Law we relate this force to the applied acceleration of the chaser, given by Eq. (15). Notice that if the applied acceleration is greater than the maximum acceleration of the chaser, the maximum acceleration is applied in the direction of the total virtual force. With the applied acceleration, the position and velocity can be obtained via integration of the equation of motion for the chaser. Total total ma a(t)= I (15) Ftotal a ota' otherwise  total  Simulation Results To illustrate the efficacy of the algorithm, we present two simple 2D simulation scenarios. As the first example, we demonstrate the differences generated by the priority weighted avoidance of the new algorithm with that of Ge and Cui,63 whose work is the baseline for the present analysis, which constrains the target and obstacles to constant velocities. In a second simulation we present a scenario more in line with the intended application, that being the Heterogeneous Expert Robots for Onorbit Servicing (HEROS) architecture,64 demonstrating a planar, onorbit rendezvous with obstacle avoidance. This simulation showcases a situation where target and obstacles have non constant velocities. Constant velocity scenario The 2D constant velocity simulation presented below is a comparison between the algorithm present in Ge and Cui63 and the algorithm derived in the present work. In the former work the repulsive potential is restricted to single obstacle interactions. The parameters (a,, a,, m, n, 7, p,) for both algorithms are equal in order to isolate the effect of the PIweighting scheme. The pertinent parameters of the simulation are as follows. The target has initial position [10 10]T and constant velocity [0.1 0.05]T. Two obstacles are in the operational space with initial positions [5 0] and [21 11]T and constant velocities [0 0.1]' and [0.05 0.065]T, respectively. The chase is of unit mass with the maximum acceleration, Amax, also equal to unity. The chase is assigned initial position [1 1]' and velocity [0.1 0]T. A FOS of 0.9 was used for the PIweighting algorithm. The results of the simulation are shown in Figure 24. 20 Target Obstl Obst2 Tatsch 15 Ge&Cui > 10 5    aI 0  j   0 5 10 15 20 25 30 X Figure 24. 2D Constant velocity simulation results. Figure 24 shows the locus of points for all entities in the operational space for the 200 seconds. Most important to note is the more direct route traversed by the chase using the P1weighting scheme when in the influence of both obstacles. As chase approaches the first obstacle in both situations, the differences in the algorithms is illustrated. The chase under Ge's and Cui's algorithm makes a more prominent maneuver with respect to the PIweighting algorithm because the former does not consider the influence of the second obstacle at all. As the chase encounters the second obstacle, the algorithms converge to the same trajectory; which is to be expected considering that in the single obstacle case, the algorithms only differ by the factor of safety. Finally we see that in the absence of obstacles the algorithms are equivalent given that they employ identical attractive potential/force interactions. Onorbit rendezvous scenario The onorbit rendezvous simulation consists of four spacecraft in Keplerian motion: target, chase, and two obstacles. The initial conditions for all the vehicles are summarized in Table 1. The objective of this scenario is to drive the chase spacecraft using the APF guidance presented above to the target spacecraft while avoiding the obstacles. Table 1. Orbit parameters for rendezvous simulation. Vehicle a [km] e i f (o Target 8000 0.1 00 00 00 Chase 9005 0.1 00 00 00 Obstacle 1 10288 0.3 00 0.07180 00 Obstacle 2 12001 0.4 00 0.71710 00 Figure 25 illustrates the result of this simulation again in the form of the locus of points for all entities in the operational space. The chase spacecraft starts radially aligned with the target at a 3.6 kilometers higher altitude. Initially the chase proceeds by thrusting in the radial direction to lower its radial velocity (and altitude), thereby increasing its orbital speed. However, it encounters the first obstacle which requires it to thrust in the negative radial direction (raising the chase's altitude) to avoid a collision. Postavoidance, the chase again thrusts as before. The second obstacle encounter requires a less drastic trajectory change. Once past all obstacles in the environment, the chase finally achieves an orbit altitude lower than that of the target and thus achieves a positive relative velocity (meaning it is catching the target without thrusting). However, the dependence on the relative position results in small accelerations still being applied to the chase until the rendezvous is achieved. 1000  S "Target 900 Obst 1 CObst 2 800 Chase 700 600 400 300 200 100 0  7160 7160 7170 7180 7190 7200 I Figure 25. Planar onorbit rendezvous results. Figure 26 is a graph of the norm of the applied acceleration of the chase spacecraft. This simulation also took into account the maximum allowable acceleration that could be achieved by the vehicle. We see that initially the chase acceleration is saturated until we encounter the first obstacle (indicated by the dip in the graph). After passing the obstacle we again see a saturation of the acceleration due to increasing lag between the chase and target spacecraft. The second dip in the graph indicates the encounter with the second obstacle. Beyond the second obstacle encounter the chase acceleration does not saturate again but rather smoothly decays to zero as it approaches the target spacecraft. 0. 0. Z 0.01 0 20 40 60 Time [s] Figure 26.. Norm of applied chaser accelerations. 80 100 Optimized Harmonic Attractive Potential Recall the two prominent flaws mentioned in the introduction: spurious local minima and suboptimal performance. The first of these is easily mitigated by utilization of harmonic potential functions, i.e., functions that satisfy Laplace's equation. By using harmonic functions for construction of the artificial potential field, the resulting potential field is guaranteed to be free of local minima. The first utilization of harmonic potential field for motion control is generally attributed to Satoh66 in the mid1980's. However the work was published in Japanese and did not receive exposure until it was reprinted in English in 1993. Masoud67 has explored both scalar and vector potential field methods using harmonic functions. For an extensive survey on these methods we refer you to the introduction of their work. More recently Waydo68 has utilized harmonic potential field methods, in the vein of potential flow analogs in hydrodynamics, for motion planning on the Cornell RoboFlag test bed. Conquering the suboptimal penalties inherent to APF methodologies without sacrificing its most attractive characteristic, that being computational simplicity, is a daunting task. However we proposed a solution using a library of harmonic function primitives (i.e., simple harmonic functions such as source, sinks, doublets in the hydrodynamic analog) to shape the force/acceleration/velocity field around the spacecraft such that it takes into account the system dynamics. A high degree of similarity can be achieved by using just two primitives, a sink vortex pair from the hydrodynamic analog, that have been optimized for sink strength a and vortex circulation F to mimic the true velocity field of a dynamic system using a genetic algorithm (GA). The reference velocity field is that for single impulse orbital intercept maneuvers (utilizing the solution of the ClohessyWiltshire equations). The parameters a and F were optimized over a 50 x 50 grid of points in terms of directional accuracy in the velocity fields. A program that generates the velocity vectors for CW transfers (with time of flight, TOF = 7/4 in this case) at all 2500 points in the grid and the velocity field for the sinkvortex pair was generated by specifying a value for a, F and subsequently quantifying the error between the two resulting vector fields. The error between the two velocity fields is quantified by the scalar error measure, E, given by Eq. (16), where i varies over all points in the grid, and v and v' are the velocities for the CW transfer and the sinkvortex field, respectively. In the ideal case (velocity fields are perfectly aligned) the angles O, are all zero, and cos(, ) equal one. Therefore, in the case of perfectly aligned vector fields, e will equal the number of points in the grid, n, in this case 2500. Using this value of E as the reference value, REF, we define the cost function, F, employed in the GA as a simple difference REF e i.e., the difference between REF and the actual value of the summation of the cosines of the angle between the velocity vectors at each point in the grid, as given by Eq. (17). The GA then minimized the value of the cost functional F for all values of a and F in their specified ranges. The genetic algorithm parameters used for the optimization are summarized inTable 2. = Zcos( ')= I (16) F = 2500 cos(O0) (17) The two parameters of a and F were restricted to range between 0.0001 and 5.0. The output of the GA for eleven runs is shown in Table 3. As is illustrated by the small standard deviations in Table 3(highlighted in blue), both parameters converged to a unique value in their respective ranges. The best fitness score achieved was 43.1029 with corresponding values of a = 0.1992 and F = 0.2941 (highlighted in red in Table 3). Figure 27 plots a comparison of the velocity fields for both the optimized sinkvortex pair (blue) and the CW field (red). Figure 27 illustrates the fidelity in the velocity field that can be achieved with just two harmonic function primitives. Table 2. Genetic algorithm statistics Parameter Name PARAMETER VALUE Population size 25 Numberof Generations 25 Probability of Mutation 0.05 Probability of Crossover 0.9 Crossover type Uniform Fitness Function Rosenbrock Number of Bits 20 Table 3. Mean and standard deviations for the two parameters in the GA optimization Run a F Cost Function Score 1 0.397 0.5309 55.2659 2 0.397 0.5309 55.2659 3 0.397 0.5309 55.2659 4 0.397 0.5309 55.2659 5 0.298 0.4451 48.8113 6 0.9852 1.3229 98.0802 7 0.1878 0.2784 43.2809 8 0.3011 0.4162 47.5251 9 0.2099 0.3171 43.5354 10 0.2838 0.3799 47.0728 11 0.1992 0.2941 43.1029 MEAN 0.3685 0.507 53.86110909 STD 0.2206 0.2881 15.52175451 Conclusions To summarize, with the results of the previous section, the guidance algorithm presented in this work is the combination of a fully dynamic obstacle avoidance algorithm coupled with an attractive potential based on harmonic functions. This form of the attractive potential is guaranteed to be local minima free. Furthermore, by using a sum of the harmonic function primitives to shape the velocity/force/acceleration field, with the intent to harness the underlying dynamics of the system, the suboptimal penalties of artificial potential field methods are mitigated. The new artificial potential functionbased guidance algorithm for dynamic environments introduced a prioritybased queuing approach to obstacle avoidance called 58 the PIweighting scheme. The PIweighting scheme gives the ability to do multiple obstacle avoidance simultaneously, and given its implicit dependence on the relative velocity between the chase (robot) and the obstacles, extends the APF guidance algorithm to general dynamic environments (nonconstant velocity targets and obstacles). *8 0.8 . 0.2 0 .. ... .. .. i ... 0.4 0.6 08 : \\\ \ 0.8 1 0.5 0 0.5 1 X Figure 27. Optimized sinkvortex pair vs. CW based vector field Through a simple 2D simulation we isolated the effects of the PIweighting scheme and demonstrated the effectiveness of the simultaneous avoidance in situations involving multiple obstacles. Furthermore, a second simulation was presented to demonstrate the efficacy of the overall APF guidance algorithm in terms of the ultimate application, the HEROS architecture. This simulation demonstrated the effectiveness of the algorithm to 0.2 ,, * t'1.4ddVGg;^ '' ',." .1 X.. .i Y Figur 27. Optimized sinkvort* ex pi vs. Cs v fr" JKSield 59 handle an onorbit rendezvous with multiple obstacles, all of which were experiencing Keplerian motion. CHAPTER 5 ARTIFICIAL POTENTIAL FUNCTION ATTITUDE GUIDANCE Thus far the spacecraft have been modeled as point masses, however at some stage we must admit the vehicles as rigid bodies (due to directionality constraints, e.g., docking axes aligment) and account for both attitude and position dynamics of the true system. Even though the attitude and position dynamics are highly coupled in the most general case, there are situations where decoupled dynamics are acceptable, and we will deal with uncoupled cases in this work, which allows independent development of an attitude guidance design. In this vein, Chapter 5 develops a novel attitude guidance algorithm utilizing the APF. The motivation for employing artificial potential functions (APF) for guidance algorithms came from the robotics literature. The literature contains volumes of position translationall) guidance algorithms, and implementations of them, but very few examples of algorithms for rigid body (6 DOF, position and attitude) guidance. Primarily the scarcity of literature on the rigid body motion planning can be attributed to two factors: the configuration space (Cspace) has dimension equal to the DOF of the robot, and configuration obstacles, configurations that result in collisions with the environment, do not have compact representation for high dimensional Cspace. Thus for all but the simplest cases, motion planning is challenging. Classically, in order to circumvent the complexities of dealing with a high dimensional Cspace, for robots with multiple DOF the attitude, or rotational, DOF have been approximated or constrained. The simplest solution for mitigating high dimensionality of the Cspace is to ignore the attitude states all together, and merely grow the robot to a sphere with radius equal to its maximum dimension. Although this does reduce the dimensionality of the Cspace, it also reduces the capabilities of the robot; for example, for inspace operations (ISO) systems, the ability to point at objects is crucial, and this approximation eliminates this possibility. An alternative employed in the literature is to hold all but one of the attitude states fixed. Again, this alleviates the issues with dimensionality, but still significantly limits the capabilities of the robot. Therefore, a new approach that allows for incorporation of all attitude states is necessary for ISO systems, preferably one that avoids calculating the Cspace and configuration obstacles. APF guidance for attitude is utilized to accomplish this task. At this point one might interject that the two deficiencies of APF methods that were discussed in the introduction to Chapter 4 would prohibit the use of the classical approach. However, in the case of attitude guidance, as opposed to position guidance, these deficiencies do not arise. The introduction of local minima which causes premature termination of the algorithm is not an issue for the attitude case because the equations of motion, as will be shown later, possess only two extrema in the state space (which are identical points in the physical space) at the equilibrium points. Suboptimal performance concerns can be dismissed by showing that the guidance law obtained from this approach belongs to a subclass of controllers obtained from a more general Lyapunov analysis. Most obviously an obstacle avoidance algorithm, dynamic or static, is no longer necessary in the attitude guidance case. In this chapter the common attitude representations and their properties are presented in order to justify the selection of the quaternions of rotation as the attitude representation scheme. Then, paralleling the procedure for developing the APF position guidance algorithm in Chapter 3, APF guidance for attitude is developed. Attitude Representation Before proceeding with the development of the APF guidance algorithms for attitude, an attitude representation scheme must be selected. For attitude representations, it is commonly known that there are no globally nonsingular 3parameter representations of the rotation group.69,70 Therefore, we look to the 4parameter sets in order to obtain global nonsingularity, and in the present work the quaternions are chosen as the attitude representation. More precisely, we are using the term quaternion as a quaternion of rotation, i.e., the norm is constrained to unity, which is precisely the column vector of the EulerRodrigues symmetric parameters.69 The primary advantage of the quaternions is that successive rotations result in successive multiplication of 4 x 4 matrices which are commutative.69,70 As a result of this computational simplicity, the quaternion attitude representation has become commonplace in computer animation71 and spacecraft contro.172'73 Furthermore, the kinematic relationship between the attitude states and the angular velocity can be written as a linear relationship,69 as given by Eq. (33). For a comprehensive survey on attitude representations see Shuster.69 The development of the attitude APF guidance is identical to the position APF guidance algorithm. However, hidden in the details is a critical difference in state variables between the two cases. For the position case, the generalized coordinates were x and the generalized velocities were x, representing the relative position and velocities in the LVLH frame, respectively. This yields the relationship between the derivative of the generalized coordinates and the generalized velocities as unity, shown by Eq. (32). For the attitude case the generalized coordinates are the elements of the quaternion i and the generalized velocities are the components of the angular velocity of the robot . The kinematic equation relating the coordinates to the velocities is no longer unity, as illustrated by Eqs. (33) and (34). Although not immediately obvious at this point what complication to the development of the guidance algorithm for attitude this presents, it will be made apparent in the subsequent development. dx =x (32) dt q= 0 q (33) 0_ = 0 3 0 0 (34) C92 O1 0 Attitude Equations of Motion The rotational dynamics of a rigid body are governed by the Euler Equations, Eq. (35) where T is the external torque and u is the vector of control torques. For the present work it will be assumed that the external torque is zero. JO +(0 x JO = Z r+u (35) The quaternion was selected to represent the attitude, which is given in terms of the unit vector in the direction of the eigenaxis of rotation, a, and angle of rotation, (p, in Eqs. (36) (38) The kinematic equations relating the angular velocity and the attitude are given by Eqs. (33) and (34) above. q =[q q4 ] (36) q= sin *a (37) q4= cos( ) (38) Again since APF guidance is based on enforcing stability about the equilibrium of the dynamical system, it is instructive to find the equilibrium solutions for the uncontrolled (u = 0) attitude equations of motion. The equilibrium states can be found by solving the following system of equations, derived from Eqs. (33) and (35): wx (J) + u = 0 (39) q x + q4 =0 (40) ) q =0 (41) From Eqs. (39) (41) we observe that the uncontrolled equilibrium state is defined as oe = 0 regardless of the steady state values of q4e and q . Continuous APF Attitude Guidance The development of continuous attitude guidance algorithms applies to situations where the actuation on board is from control moment gyros (CMG) clusters, reaction wheels, or similar continuous input devices. The control is derived utilizing the general procedure for APF methods discussed in Chapter 3. This general procedure will be repeated here for completeness. Choose a scalar function V(x) of the states x, such that V is positive definite for all x # x, where x, is the equilibrium state. Furthermore, we constrain the control u such that its action is always in the direction directly opposite the gradient of V(x). Employing this control ensures the total derivative of V is negative definite for all x xe, thereby we are guaranteed by Lyapunov's 2nd theorem that x > _x as t > o. Following the above procedure, consider the positive definite scalar potential function V(q,a )= ( qTPq + Co'Q), where q is the vector component of the quaternion and P and Q are positive definite and diagonal (for simplicity) matrices. The utilization of the vector component of the quaternion only is justified by the existence of the unity norm constraint, by which the scalar part is computable from knowledge of the vector component. Implementing a continuous control by constraining the attitude rates, under the action of the control (q ) to be proportional and directly opposite to the resultant of the gradient of V with respect to q and c, this can be thought of as a pseudosteepest decent approach, given by o+ = k(VV +VVV) (42) where V V = Pq and VmV = QQ. In order to obtain the control input Aq necessary to achieve this condition, bring the attitude rates before control application, q to the right hand side of Eq. (42) resulting in Eq. (43). Since the attitude rates q are noncausal, the control must be implemented via the angular velocity. Using the kinematic equation to relate 4_ and 0o and assuming the control input to the attitude rates is equivalent to the control input in the angular velocity states, i.e., Aq4 = IAo, yields the control u in Eq. (71), where G, is given by Eq. (45). Aq = kPq kQm q (43) u = kPq G, (44) G, = kQ+( q41 ) (45) 2 It is instructive at this point to address the issue of suboptimality for APF methods. The control defined in Eq. (45) is a special case of a more general quaternion feedback regulator, shown in Eq. (46), developed in a seminal paper by Bong Wie,75 et al. This larger class of quaternion feedback regulators was shown, with proper gain selection, to provide neareigenaxis (nearoptimal) rotations with guaranteed global stability. For a precise treatment of the general problem we refer you to Ref. 75. We will prove the stability of the special case in the next section. u= (nJ+D)mKq (46) Before proceeding with the stability analysis, it is instructive at this point to evaluate the closedloop equilibrium states now that we have defined a control law. Returning to the equilibrium equations in Eqs.(39) (41), we find that for u = 0 to be satisfied at the closedloop equilibrium state, both the angular velocity and the vector part of the quaternion must be zero (i.e., e, = qe = 0). Therefore, the closedloop system has two equilibrium states as shown below; however, they represent the same state in the physical space. e = 0 (47) q (48) 41 (49) Stability Analysis The following theorems, along with conditions set forth in the subsequent lemmas, establish the global asymptotic stability of the equilibrium state, as defined previously. THEOREM 1: The closedloop system defined by Eqs.(35) and (44) is globally asymptotically stable if M is positive definite, where M is defined by Eq. (50). M = (50) ( (Q +1) P k2 PROOF: Consider the candidate Lyapunov function V(q, c) = ( Pq + or Jo), where J is the inertia dyadic of the rigid body. It is important to note that the inertia matrix is used in the place of Q due to an elegant simplification that results during the stability proof. The expression for the total derivative of V is shown in Eq. (51). Substituting for Jo from Eq. (35), with the external torque r = 0, and for q from the control constraint given by Eq. (42) we obtain Eq. (52). Substituting the control given by Eq. (44) and simplifying yields Eq. (53) (NOTE: the simplification resulting from the use of J as the derivative matrix gain is the elimination of the cross product term in Eq. (52). V = qT pq + CJor j6( (51) = q (kPqkQ)+T ( Jwo+u) (52) = kq P2q kq PQ) OkPq o Gv (53 (53) Since V is a scalar function, we may transpose any term in the above expression while maintaining its value. Therefore, transposing the second term in Eq. (53) and 68 combining with the third term leads to the expression in Eq. (54). Rearranging yields a compact expression for the derivative in a quadratic form given by Eq. (55). S= kq P2q T (kQP + P)q_ G(54) k q k (_Q + 1) p _ qP 10)] kP2 k 2 (55) Now, observe that the matrix in Eq. (55) is precisely M. Therefore, if M is shown to be positive definite, then V is negative definite and by Lyapunov's 2nd theorem global asymptotic stability is guaranteed. Note, Gv being positive definite is a necessary condition for the positive definiteness of M. LEMMA 1: The derivative gain matrix, Gv, is positive definite if and only if 1 kQ, > . 2 PROOF: The derivative gain matrix is precisely given by Eq. (56) below. Positive definiteness requires that the principle minors and the determinant are positive. The ith principle minor of Gv is shown in Eq. (57) with i = 1, 2,3 and i # j # k. The principle minors are positive ifEq. (58) holds. Given that 1 < q4 < 1, Eq. (58) implies that the 1 necessary condition for the positive definiteness of Gv is kQ, > . 2 kQ + 4 q3 q2 kQ _+ 2 2 2 G 1= kQ2 + 4 2 2 2 q2 q__ 2kQ3 + 2 2 2 (56) (56) PM, = kQ +4 kQk + q(5)+ : 2 2 2 + ) LI (57) kQ, + q4>0 2 (58) The determinant of G,, given by Eq. (59), is always positive for positive diagonal 1 elements of G,. Therefore, kQ, > is a necessary and sufficient condition for the 2 positive definiteness of G,. 3 3 det (G,) = Iv,, +4I Gv,, (59) 1 =1 y = 1 2 It is instructive to point out that, in general, that all diagonal elements kQ +4 2 could be negative and still satisfy the positivity of the principal minors as required by Eq. (57). Since Q was assume diagonal and positive definite, this implies that k < 0. It will be shown that k > 0 is required for M to be positive definite, and therefore, the condition 1 kQ, > will be used. 2 LEMMA 2: A block matrix S, given by Eq. (60) is positive definite if and only if the matrices A and R = D BA B are positive definite. S =[ A] (60) PROOF: Start by executing block elementary column operation on S in order to make the resulting matrix, T shown in Eq. (61), block upper triangular. From Eq. (61) we see that the upper left block entry is exactly R Also, from its upper triangular form, T is positive definite if and only if A and R = D BA 1B are positive definite. Since elementary operations, similarity transformations in general, do not alter the eigenvalues, it is implied that if T is positive definite only if S is positive definite and vice versa. T = [DBA'B B T= A ] (61) 0 A 1 THEOREM2: For kQ, > , M is positive definite. PROOF: In Lemma 2 let S = M (Note: M has the necessary structure) which gives the matrices A, B, and D from Eq. (60) as A = kP2 B= (Q+1)P (62) 2 D=G, Using the result of Lemma 2, if A and R = D BA 'B are positive definite then M is positive definite. In this case A = kP2 which needs to be positive definite, which implies that k > 0. Furthermore, since P is a diagonal matrix so is A, which implies that the inverse of A is just the reciprocal of the diagonal elements A' =diag (Y,' A 1A. Y This gives R the form shown in Eq. (63). B2 D A D12 D13 R= D12 D 2 D23 (63) A2 D3 D23 A A 3 R has the same structure as GV, hence it follows from Lemma I that R will be positive definite if the condition D I > 0 is satisfied. Assuming that A > 1 (i.e., p2 > /lk) and substituting for A, B, and D the condition for positive definiteness of R becomes the expression in Eq. (64). kQ + 4 k2( +1)22 > 0 (64) Let the left hand side of the inequality in Eq. (64) be represented by W(P,) = c a2 as shown in Figure 1, where 2c = 2kQ, +q4 and 4a = k2(Q, +1)2. Therefore, R is positive definite for all values of ], such that W is positive. W is the equation of a parabola, with c characterizing the W intercept and a governing the width. Examining the range of the two constants a and c, we find that a > 0 always, which implies W is always opening downward, and c > 0 from kQ, >1/2 and 1< q4 < 1. This implies that for all choices of kQ, > 1/2 a portion of W lies above the ]4 axis as can be seen in Figure 1 below. Therefore, for kQ > 1/2 we can find a range of P4 such that R is positive definite and therefore M is positive definite. Furthermore, the range on ], is given by previously, meaning ], must be selected from the positive side of this range. It is important to note that k is a free parameter in the sense that it is not explicitly specified in the selection of the control gains. Consequently, the gains may be chosen within their range of stability specified by lemmas and theorems and a k found after the fact to guarantee the stability. The restriction that 0 < I < c la along with the other constraints imposed on the controller parameters from Theorems 1 and 2 and Lemmas 1 and 2 make selection of the control gains Gp and GD cumbersome. To circumvent this difficulty we restrict our control gains to a region of interest and employ an optimization algorithm, specifically a genetic algorithm discussed in the next section, to determine the control gains G, and GD. W (c/a) 1 2 Figure 28. Inequality of Eq. (64) +(c/a)'2 Pi stable Pi Stability of Dynamical System The convergence of the vector component of the quaternion and the angular velocity to zero for the controller derived in the present work is guaranteed from the analysis of the previous section. One would ask at this point as to the behavior of q4 in the dynamical system. For a general perspective on the nature of the equilibrium states given by Eqs. (39) (41), look at the linearization of the closed loop dynamical system, 73 given by Eqs. (35) and (44), around both equilibrium states. The behavior of the dynamical system in the neighborhood of these states is elucidated by the eigenvalues of these linearizations. It can be shown that for q4 = +1 the linearization is stable, implying that the spectrum of the eigenvalues, o(A) < 0. For q4 = 1 the dynamical system is unstable, meaning there exist A > 0. Therefore, regardless of the initial value of q4, it will converge to q4 = +1, which depending on the application may be viewed as a blemish on the efficiency of the controller. A plot illustrating the behavior of q4 when the initial states are very close to desired (equilibrium) state in the physical space, but the unstable equilibrium in the state space is shown below. 1 0.8  0.6 0.2  0.2 S20 Time [s 5 10 15 20 25 Time [s] Figure 29. Behavior of q4 Optimization To aid in the selection of the control gains G, and GD a genetic algorithm was used to optimize the gains with respect to a cost function, F given by Eq. (17) that is a function of the error in the states X [= ) ( ] and the control effort u. The matrices Hx and HD are weight functions for the state error and control states respectively, and for the optimization performed here both were assigned the appropriated dimensioned identity matrix, meaning equal weighting was used for the control and state errors. Since the objective of the control is to bring the state vector to the equilibrium state defined by Eq. (75) (76), the error states are equivalent to states themselves [r T [r Tr F = XerrHX +ul u (65) 2  2  The genetic algorithm parameters used for the optimization are summarized in Table 4. Table 4. Genetic algorithm statistics Parameter Name PARAMETER VALUE Population size 25 Numberof Generations 25 Probability of Mutation 0.05 Probability of Crossover 0.9 Crossover type Uniform Fitness Function Rosenbrock Number of Bits 20 The fitness function evaluated the contribution of the error in the states, q, q4 and ca and the control effort over the first 50 seconds of the simulation given by the numerical example in the next section. The control matrices G, and GD are 3x3 diagonal matrices, and therefore six parameters were optimized. The mean and standard deviation for these six parameters for ten runs of the genetic algorithm are summarized in Table 5. Table 5. Mean and standard deviation for the six parameters in the GA optimization Parameter Mean Std. Deviation Gpi 1.1083 0.1083 Gp2 1.0111 0.0101 Gp3 1.1703 0.1293 GD1 0.5253 0.0218 GD2 0.5527 0.0433 GD3 0.8038 0.1094 Cost Function 49.5967 0.4197 The three parameters of G, were restricted to range between 1.0 and 2.0, and those of GD were restricted to 0.5 to 1.0. As is illustrated by the small standard deviations in Table 5, in general all parameters converged to a unique value in their respective ranges. Numerical Example In this section, we use the numerical example given in Wen and KreutzDelgado70 to compare our controller to their classical PD controller. The objective of the numerical example is to maneuver a spacecraft, represented by its inertia matrix, J, from an initial orientation and angular velocity to the equilibrium state given by Eqs. (75) (76),, with q4 = +1. The controller in Ref. 70 is defined as shown in Eq. (66), where kp = 4 and k =8. CO) 0..4 0 5 10 15 20 Time (s) Figure 30. Time response of vector part of quaternion From Ref. 70, the spacecraft's inertia matrix is given by by Eq. (77) and the initial orientation in axisangle representation is a = [0.4896 0.2052 0.8480]T and 0 = 2.4648 radians. The spacecraft's initial angular velocity is = [0.1 0.3 0.5] The gains for the controller developed in the present work are listed in Eqs. (68) (69). u = k qk, (66) 1.0 0 0 J= 0 0.63 0 (67) 0 0 0.85 1.1083 GP = 0 0 0.5253 GD= 0 0 0 1.0111 0 0 0.5527 0 0 0 1.1703 0 0 0.8038 (68) (69) Figure 30 illustrates the response of the three components of the vector part of the quaternion, q, and Figure 31 shows the response of the scalar part of the quaternion, q4. Figure 32 and Figure 33 show the response of the three components of the angular velocity and the control input, respectively. As demonstrated by these plots, the controller developed in this paper is capable of performing the stated maneuver. In fact, the maneuver was accomplished at a significant reduction in control effort. It is important to note however, that the gains were arbitrarily selected in Ref. 70 and no optimization was investigated. 5 10 15 Time (s) Figure 31. Time response of scalar part of quaternion 20 25 78 0.6  w 1wen ...... W2wen 0 .4    .2n w W 3wen o W1 o 0 .2    02 . . .. 0.4   0 0.4  0 0.8 0 5 10 15 20 25 Time (s) Figure 32. Time response of angular velocity The effective performance of the controller was demonstrated through comparison with a standard PD controller by a numerical example. It is important to note that even though the gain matrices components were in general smaller that the gains used in the PD controller simulation, the performance achieved was superior, as is illustrated by Figure 30 through Figure 32. Furthermore, the control effort to achieve this improved performance was considerably less, as can be seen in Figure 33. It should be stated here that the convergence of q4 to +1 introduces a possible inefficiency for attitude maneuvers. If initially q4 is closer to 1 than +1, then it would be more efficient to drive the attitude toward the equilibrium state at q4 = 1. However, our controller drives the attitude to the equilibrium state with q4 = +1, regardless of the initial value of the scalar part of the quaternion. One possibility to mitigate this inefficiency is to make the gain matrix G, in Eq. (71) proportional to the signum(q4), as done in Wie, et. al.73 0 5 10 15 Time [sl Figure 33. Time response of control effort 20 25 Impulsive APF Attitude Guidance Impulsive attitude guidance is applicable in those situations where the actuation is provided by thrusters as opposed to the devices used in the continuous case, which with the growing sophistication and maturity of microNewton capable microthrusters is becoming a popular option for attitude control systems (ACS) on board micro and nano class satellites. 80 Paralleling the development in previous section, let the artificial potential function for this case be V(q) = q Qq where q is the vector component of the quaterion q. Note, it is acceptable to use the vector portion of i only due to the constraint on quaternions of rotation, i.e., 11711 = 1. The total derivative of the potential function is given by V = (VqV) q, and the control is defined by Eq. (71) where is such that the rate of the vector part of the quaternion after a control impulse q is directly opposite to the local gradient of the potential qo = kVqV Defining the control in this manner ensures the derivative of the potential after a control impulse is negative definite, as given by Eq. (70), and therefore, by Lyapunov's 2nd method q is guaranteed to converge to its equilibrium value. V'(x)= k(VV) (VV) < 0 (70) o, r'<0 Control = (71) \8, V> > 0 S _= 0 0 =kVV (72) At this point, the difference in the relationship between the generalized coordinates and velocities become apparent. Since q are not causal states, the control must be applied to the angular velocity states. Utilizing the kinematic relationship given by Eq. (33) the control constraint can be transformed to a form in terms of the angular velocity 0, given by Eq. (72), while maintaining the constraint on q Substituting the gradient Vq = Qq and premultiplying both sides of Eq. (72) by the inverse of the coefficient matrix of 0,, yields Eq. (73). Thus the control is applied to the angular velocity states such that an impulsive angular velocity correction Ao = 0, g, that constrains the rate of the vector component of the quaternion to be tangent to the local gradient of the potential. o0 =k J Qq (73) q 0  Another difference for APF guidance development between the position and attitude cases is the proof of the stability of the equilibrium states. The equilibrium was the origin in the position case, with the equilibrium derived from the orbital relative motion equations of motion. For the attitude we must return to Eqs. (39) (41) now that the control has been defined in order to obtain the controlled (u_ 0) equilibria. Just as before, the equilibrium equations show that for u = 0 to be satisfied at the equilibrium solutions both the angular velocity and the vector part of the quaternion must be zero (i.e., o)e = q = 0 ). Again there are two equilibrium states (both the same point in the physical space) for the controlled system, restated in Eqs. (74) (76). Since no differentiation between the equilibrium solutions is necessary in the convergence analysis above, the states will converge to the equilibrium solution to which it is initially closest. co =0 (74) q = 0 (75) q4. = +1 (76) Numerical Example: Impulsebased APF Attitude Guidance To illustrate the efficacy of the impulse based APF guidance algorithm for attitude, a numerical example is considered. A satellite in orbit has an initial orientation given in axis/angle representation as a = [0.4896 0.2032 0.8480]T and 0 = 2.4648 radians, and its initial angular velocity is given by = [0.5 0.8 0.9]' radians per second. The central principle inertia matrix representing the spacecraft is given by Eq. (77). The APF parameters used were k = 0.1 and Q = I, where I is the 3 x 3 identity matrix. The objective of the APF guidance algorithm is to bring the spacecraft with its initial parameters to the equilibrium state [qT q4 'T T = 0 +1 01. 1.0 0 0 J= 0.63 0 (77) 0 0 0.85 0.7 . 0.6 0.1 S 0 A.1 1 05 10 13 20 25 Time Is] Figure 34. Components of the vector part of the quaternion vs. time .0.2 0.  0. 0.8 0 5 10 15 20 25 Time [s] Figure 35. Components of the angular velocity vs. time 0.9 as I 8 7  I R 0.5 6 0.4 0 5 10 s1 20 25 Time [s] Figure 36. Scalar part of the quaternion vs. time 84 Figure 34 and Figure 35 show q and a versus time, respectively. Figure 36 illustrates the trajectory followed by q4 as a function of time. The maneuver duration was 25 seconds and a total impulsive angular velocity change Aca,,, = 1.4034 radian per second was necessary to bring the <1 10 3. The impulsive angular velocity inputs are illustrated by Figure 37, and it can be observed that only three impulses totaling 1.4034 radians per second were necessary for the maneuver. Solely for the purpose of illustrating the convergence to the equilibrium point initially closer, the numerical example was repeated with identical parameters, except the initial value of q4 = 4.4648 radians. Figure 38 depicts the trajectory followed by the scalar part of the quaternion in this case, and it is observed that indeed it converges to 1. O 5 10 15 20 25 0 5 10 15 20 25 4, 0 ,T[, Time FR1 O". TI__,. .2 II Time [s] Figure 37. Angular velocity impulses vs. time for the numerical example An important caveat to the results presented here the relatively fast convergence of the attitude and angular velocities are direct consequences of a small inertia matrix in Eq. (77) and small initial angular velocities. However, the results here are for illustration of the efficacy of the algorithms to bring the states to their equilibrium values. A more precise analysis of the performance of these algorithms is necessary utilizing more practical inertias and initial parameters. 0.65 0,7 0.75 0.8 D.85 0.9 0.95 Time [s] Figure 38. Scalar part of the quaternion vs. time for minus one convergence case Conclusions In this chapter two original attitude guidance algorithms were derived from an artificial potential function approach. In the continuous guidance case, the algorithm was shown to belong to a more general class of quaternion feedback regulator to eliminate concerns about suboptimal performance inherent in APF methods, and the algorithm was proven to be globally asymptotically stable. The controller parameters were obtained using a genetic algorithm in order to mitigate the difficulty in gain selection and to obtain a degree of optimization with respect to state error and control effort. The GA obtained convergence for all six parameters in their respective ranges, as was illustrated by the small standard deviation of each parameter. Furthermore, an impulsive guidance approach was developed and shown to be effective for spacecraft with thruster based ACS. With these algorithms in place, a more complete development of APF guidance that synthesized both position and attitude APF algorithms for rigid bodies is pertinent. These are the topics under current investigation and will be the directions pursued in the future. In the next chapter, we conclude this dissertation with some recommendations for future work and our concluding remarks. CHAPTER 6 CONCLUSION AND FUTURE WORK Conclusion The current direction of the space industry is toward utilization of autonomous space robotic systems for inspace operations (ISO) and planetary surface exploration (PSE). The state of the art outlined in Chapter 1 illustrates few achievements in real world implementations and test bed environments for ISO space robotics. As a result research in autonomous space robotic systems is both timely and necessitated. With this as the stimulus we explore autonomous guidance techniques for advancing the state of the art ISO space robotic systems. As a first attempt the classical orbit mechanics methods were explored in Chapter 2. These methods were shown to be insufficient, primarily due to the inability to incorporate obstacle avoidance, for a large class of ISO mission. In Chapter 3 we turned to the robotics literature for alternative approaches and discovered artificial potential function (APF) guidance to be a feasible solution, provided the deficiencies inherent in these methods are overcome; the lack of dynamic obstacle avoidance, spurious local minima, and suboptimal performance. Chapter 4 developed the key elements that allow APF guidance to be an effective solution for ISO missions. A novel dynamic obstacle avoidance algorithm is generated that is capable of handling general dynamic environments. Furthermore, through the use of a sum of harmonic function primitives we can guarantee no introduction of local minima and mitigate the suboptimal penalties of APF methods. At some point the point mass assumption for the spacecraft becomes unacceptable and instead must be admitted as rigid bodies. As a result, consideration of the attitude in introduced. In Chapter 5 the classical APF framework is utilized to develop two novel attitude guidance algorithms, one for continuous and one for impulsive type actuation. From characteristics of the attitude control problem and proper potential function selection eliminates the inherent deficiencies of an APF approach. Both controller were shown to perform effectively for their appropriate application. Future Work The solution method presented in this dissertation is an effective solution for autonomous guidance of inspace operations mission, such as assembly, inspection, and servicing. However, further development of this approach is possible in several areas, as outlined below. Harmonic Functions o Utilization of more complex primitives o Use additional primitives in sum o Eliminate the use of offline optimization so algorithm could be used for real time guidance Reference Velocity Field o Generate full nonlinear velocity field (Lambert solution field) o Include perturbing effects in reference field Drag J2 Solar radiation 89 System Dynamics o Fuse attitude and position guidance algorithms for couple attitude/position dynamics Research in these areas is the focus of ongoing post graduate study. 