University of Southern California Dissertations and Theses

Optimal guidance trajectories for proximity maneuvering and close approach with a tumbling resident space object under high fidelity J₂ and quadratic drag perturbation model

(USC Thesis Other)

Optimal guidance trajectories for proximity maneuvering and close approach with a tumbling resident space object under high fidelity J₂ and quadratic drag perturbation model

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content Optimal Guidance Trajectories for Proximity Maneuvering and Close Approach with a Tumbling Resident Space Object under High Fidelity J 2 and Quadratic Drag Perturbation Model by Parv Patel A Dissertation Presented to the FACULTY OF THE GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulllment of the Requirements for the Degree DOCTOR OF PHILOSOPHY (Astronautical Engineering) December 2021 Copyright 2021 Parv Patel Dedication To Mom and Dad, The reason of what I am today. Thank you for your great support and continuous care. To my sisters and brothers, I am really grateful for all your love and encouragement. To Science, Thank you for enlightening me. ii Acknowledgements Every Ph.D. is the result of not only a single person's eorts, but the broader belief and striving of an entire civilization. This one in particular is the product of many things, so I'll try to keep it brief. First and foremost, I owe an enormous thank you to my advisor and my long term mentor, Dr. Bogdan Udrea, for taking me on as a student, nurturing my interest in arcane and varied topics, and pushing me to achieve at a level I had never thought possible. His guidance, teaching and encouragement carried me throughout graduate school, and will continue to carry me throughout my career. I would like to thank my parents, Keyur and Dharmita, whose enthusiasm, support, and belief in my abilities kept me going in rough times. My sister, Pratha and all my cousins, Vishva, Hayaa, Katha and Srusty, their extensive support throughout my time as a student and for always being their for all my ups and downs. To my second family, thank you Jigoomama, Chiragmama, Dharmimami and Vrundamami for all your love and unfailing emotional support. There are many other people in the USC community who were tremendously helpful in the production of this dissertation, too many to list. I would like to specically thank each of the members of my PhD committee for volunteering their time and knowledge to aid in the writing of this dissertation. Dr. Dan Erwin for always guiding me during my toughest times and lastly Dell Cuason for always taking care of me. To all my friends who have been an integral part of my growth, I am very thankful to have learned so much from each and every one of you. Thank you to Ankush Sahni, Navroz Lamba, iii Aman Sharma, Rizwan Patel and Maggie Schlossberg for always encouraging me to pursue the unknown. A special thanks to Peter Zhang for showing me what true work ethic is, there is so much more I have yet to learn from you. And nally my lovely and beautiful niece for always putting a smile on my face. This is for you Tyaagi! iv Table of Contents Dedication ii Acknowledgements iii List Of Figures vii Abstract xi Chapter 1: Introduction 1 1.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Motivations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.3 Literature Review . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.4 Organization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 Chapter 2: Relative Motion Models 8 2.1 Coordinate Frames . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.1.1 ECI - Earth Centered Inertial Coordinate Frame . . . . . . . . . . . . . . . 9 2.1.2 ^ r ^ ^ i - Spherical Coordinate Frame . . . . . . . . . . . . . . . . . . . . . 9 2.1.3 LVLH - Local Vertical Local Horizontal Coordinate Frame . . . . . . . . . . 10 2.1.4 ^ x ^ y ^ z - Curvilinear Coordinate Frame . . . . . . . . . . . . . . . . . . . 11 2.2 Clohessy-Wiltshire Relative Motion Model . . . . . . . . . . . . . . . . . . . . . . . 12 2.3 High Fidelity J 2 Relative Motion Model . . . . . . . . . . . . . . . . . . . . . . . . 17 2.3.1 Time AveragedrJ 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 2.3.2 Correcting the orbital period of the Reference orbit . . . . . . . . . . . . . 21 2.3.3 Correcting the Reference orbit for nodal drift . . . . . . . . . . . . . . . . . 22 2.3.4 Time Averaged J 2 motion model with corrections . . . . . . . . . . . . . . 23 2.4 Quadratic Drag Relative Motion Model . . . . . . . . . . . . . . . . . . . . . . . . 24 2.5 Full Linear Model with Corrected J 2 and Quadratic Drag Dynamics . . . . . . . . 30 Chapter 3: Eects Of Perturbation 35 3.1 Eects of the J 2 Disturbance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 3.2 Eects of Drag Disturbance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 3.3 Validating the Full Linear Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 3.4 Comparisons with Nonlinear Model . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 Chapter 4: Control Problem Formulation 68 4.0.1 Optimization Problem Setup . . . . . . . . . . . . . . . . . . . . . . . . . . 70 4.1 Minimum Control Solution for Full Linear Model . . . . . . . . . . . . . . . . . . . 72 4.1.1 Eect of Non Corrected J 2 Perturbation . . . . . . . . . . . . . . . . . . . . 75 4.1.2 Eect of Inclination on the Perturbed Optimal Trajectories . . . . . . . . . 76 v 4.2 Minimum Control Solution for Hill's EOMs . . . . . . . . . . . . . . . . . . . . . . 77 4.3 Cost Function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 4.4 Eect of Altitude and Maximum Allowable Thrust . . . . . . . . . . . . . . . . . . 84 Chapter 5: Linear Quadratic Regulator 87 5.0.1 Choosing Q and R matrices . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 5.0.2 Numerical Solutions to the LQR problem . . . . . . . . . . . . . . . . . . . 90 5.1 LQR Control Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 5.1.1 Variable t f . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 5.1.2 Comparing v . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 Chapter 6: Neuromorphic Controller 97 6.0.1 Neural Engineering Framework . . . . . . . . . . . . . . . . . . . . . . . . . 102 6.1 Principle 1 { Representation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 6.2 Principle 2 { Transformation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 6.3 Principle 3 { Dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108 Chapter 7: Conclusions 112 7.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112 7.2 Contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113 7.3 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 References 120 vi List Of Figures 2.1 ECI Coordinate Frame and ^ r ^ ^ i Spherical Body-Fixed Coordinate Frame . . . 9 2.2 LVLH Body-Fixed Coordinate Frame . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.3 ^ x ^ y ^ z Curvilinear Body-Fixed Coordinate Frame . . . . . . . . . . . . . . . . . 11 2.4 x - relative position between two spacecrafts dened in LVLH . . . . . . . . . . . . 12 3.1 Relative motion trajectories with i ref = 70 for 2-orbital period . . . . . . . . . . . 37 3.2 Relative position errors between Hill's Model and Time-averaged J 2 with correc- tions Model (i ref = 70 and 2-orbital period) . . . . . . . . . . . . . . . . . . . . . 38 3.3 Relative position errors between Time-averaged J 2 with corrections Model and J 2 without corrections [Eq. 2.39] (i ref = 70 and 2-orbital period) . . . . . . . . . . . 39 3.4 Relative motion trajectories with i ref = 0 for 2-orbital period . . . . . . . . . . . . 40 3.5 Relative position errors between Hill's Model and Time-averaged J 2 with correc- tions Model (i ref = 0 and 2-orbital period) . . . . . . . . . . . . . . . . . . . . . . 41 3.6 Relative position errors between Time-averaged J 2 with corrections Model and J 2 without corrections [Eq. 2.39] (i ref = 0 and 2-orbital period) . . . . . . . . . . . . 42 3.7 Relative motion trajectories at 900-km altitude with i ref = 70 for 2-orbital period 44 3.8 Relative position errors between Quadratic Drag Model and Hill's Model (i ref = 70 with 900-km Altitude for 2-orbital period) . . . . . . . . . . . . . . . . . . . . 45 3.9 Relative motion trajectories -i ref = 70 with I.C. = [1; 0;1] for 2-orbital period (900-km altitude) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 3.10 ^ x ^ z view for I.C. = [1; 0;1] for 2-orbital period (900-km altitude) - Drift view 46 3.11 Relative position errors between Quadratic Drag Model and Hill's Model for I.C. = [1; 0;1] (i ref = 70 with 900-km Altitude for 2-orbital period) . . . . . . . . . 47 vii 3.12 Relative motion trajectories with Quadratic Drag Model at multiple altitudes for I.C. = [1; 0;1] (i ref = 70 for 2-orbital period) . . . . . . . . . . . . . . . . . . . 48 3.13 Zoomed view of Relative motion trajectories with Quadratic Drag Model at multi- ple altitudes for I.C. = [1; 0;1] (i ref = 70 for 2-orbital period) . . . . . . . . . 49 3.14 Relative position errors with Quadratic Drag Model for dierent altitudes for I.C. = [1; 0;1] (i ref = 70 for 2-orbital period) . . . . . . . . . . . . . . . . . . . . . 50 3.15 Relative motion trajectories -i ref = 70 with I.C. = [1; 0;1] for 7-orbital period (900-km altitude) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 3.16 ^ x ^ z view for I.C. = [1; 0;1] for 7-orbital period (900-km altitude) - Drift view 52 3.17 Relative position errors between Quadratic Drag Model and Hill's Model for I.C. = [1; 0;1] (i ref = 70 with 900-km Altitude for 7-orbital period) . . . . . . . . . 52 3.18 Relative motion trajectories for I.C. = [1; 0; 0] (i ref = 70 with 900-km Altitude for 2-orbital period) and = = 0 . . . . . . . . . . . . . . . . . . . . . . . . . . 53 3.19 Relative position errors between Hill's Model and Quadratic Drag Model for I.C. = [1; 0; 0] (i ref = 70 with 900-km Altitude for 2-orbital period) and = = 0 54 3.20 Relative position errors between Hill's Model and Quadratic Drag Model for I.C. = [1; 0; 0] (i ref = 70 for 2-orbital period) with = = 0 and h a =h p = 902-km 55 3.21 Relative motion trajectories for Full Linear Model with both J 2 and quadratic drag perturbation for 2 orbital-period (i ref = 70 and 900-km altitude) with IC = [1; 0; 0] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 3.22 Relative position errors between Hill's Model and Full Linear Model with both J 2 and quadratic drag perturbation for 2 orbital-period (i ref = 70 and 900-km altitude) with IC = [1; 0; 0] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 3.23 Relative motion trajectories for Full Linear Model with both J 2 and quadratic drag perturbation for 2 orbital-period (i ref = 0 and 900-km altitude) with IC = [1; 0; 0] and = = 0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 3.24 Relative position errors between Hill's Model and Full Linear Model with both J 2 and quadratic drag perturbation for 2 orbital-period (i ref = 0 and 900-km altitude) with IC = [1; 0; 0] and = = 0 . . . . . . . . . . . . . . . . . . . . . . 59 3.25 Relative motion trajectories for Full Linear Model with both J 2 and quadratic drag perturbation for 2 orbital-period (i ref = 0 and 900-km altitude) with IC = [1; 0; 0] compared to Nonlinear Dynamics Model . . . . . . . . . . . . . . . . . . . 60 3.26 Relative position errors between Nonlinear Dynamics Model and Full Linear Model with both J 2 and quadratic drag perturbation for 2 orbital-period (i ref = 0 and 900-km altitude) with IC = [1; 0; 0] . . . . . . . . . . . . . . . . . . . . . . . . . . 61 viii 3.27 Relative motion trajectories for Full Linear Model with both J 2 and quadratic drag perturbation for 2 orbital-period (i ref = 70 and 900-km altitude) with IC = [1; 0; 0] compared to Nonlinear Dynamics Model . . . . . . . . . . . . . . . . . . . 62 3.28 Relative position errors between Nonlinear Dynamics Model and Full Linear Model with both J 2 and quadratic drag perturbation for 2 orbital-period (i ref = 70 and 900-km altitude) with IC = [1; 0; 0] . . . . . . . . . . . . . . . . . . . . . . . . . . 63 3.29 Relative motion trajectories for Full Linear Model with both J 2 and quadratic drag perturbation for 5 orbital-period (i ref = 0 and 900-km altitude) with IC = [1; 0; 0] compared to Nonlinear Dynamics Model . . . . . . . . . . . . . . . . . . . 64 3.30 Relative position errors between Nonlinear Dynamics Model and Full Linear Model with both J 2 and quadratic drag perturbation for 5 orbital-period (i ref = 0 and 900-km altitude) with IC = [1; 0; 0] . . . . . . . . . . . . . . . . . . . . . . . . . . 65 3.31 Relative motion trajectories for Full Linear Model with both J 2 and quadratic drag perturbation for 5 orbital-period (i ref = 70 and 900-km altitude) with IC = [1; 0; 0] compared to Nonlinear Dynamics Model . . . . . . . . . . . . . . . . . . . 66 3.32 Relative position errors between Nonlinear Dynamics Model and Full Linear Model with both J 2 and quadratic drag perturbation for 5 orbital-period (i ref = 70 and 900-km altitude) with IC = [1; 0; 0] . . . . . . . . . . . . . . . . . . . . . . . . . . 67 4.1 Space vehicle thruster conguration as l 1 norm . . . . . . . . . . . . . . . . . . . . 69 4.2 Optimal trajectory for the Full Linear Model (i ref = 70 ) . . . . . . . . . . . . . . 72 4.3 Planar view of the optimal trajectory for Full Linear Model (i ref = 70 ) . . . . . 73 4.4 Position and Velocity proles for the Full Linear Model (i ref = 70 ) . . . . . . . . 74 4.5 Minimum propellant control proles for the Full Linear Model (i ref = 70 ) . . . . 74 4.6 Optimal guidance trajectory for dierent models (i ref = 70 ) . . . . . . . . . . . . 76 4.7 Optimal trajectory for the full linear model for dierent inclination . . . . . . . . . 76 4.8 Optimal trajectory for the Hill's Equation and the Full Linear Model (i ref = 0 and = = 0). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 4.9 Minimum propellant control proles for the Full Linear Model for i ref = 0 and = = 0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 4.10 Minimum propellant control proles for the Hill's Equation . . . . . . . . . . . . . 78 4.11 Optimal trajectory for the Full Linear Model (i ref = 70 ) . . . . . . . . . . . . . . 80 4.12 Minimum propellant control proles for the Full Linear Model for i ref = 70 with 6 eective controllers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 ix 4.13 Minimum propellant control proles for the Full Linear Model for i ref = 70 with 3 eective controllers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 4.14 Optimal trajectory for the Hill's Relative Motion Model (i ref = 0 ) . . . . . . . . . 82 4.15 Minimum propellant control proles for the Hill's Model for i ref = 0 with 6 eective controllers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 4.16 Minimum propellant control proles for the Hill's Model for i ref = 0 with 3 eective controllers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 4.17 Minimum propellant trajectories for variable altitude and dierent time-to-destination constraint . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 4.18 Minimum propellant trajectories for variable altitude and dierent time-to-destination constraint . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 5.1 Schematic of the liner quadratic regulator (LQR) for optimal full-state feedback. The optimal controller for a linear system given measurements of the full state, y = x is given by proportional control u =K r x where K r is optimal gain at a given time t . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 5.2 LQR Optimal trajectory for CHW Model (1 orbital period) . . . . . . . . . . . . . 93 5.3 Position and Velocity proles for CHW Model (1 orbital period) . . . . . . . . . . 93 5.4 LQR optimal trajectory control prole for CHW Model (1 orbital period) . . . . . 94 5.5 LQR Optimal trajectory for CHW Model (2.572 orbital period) . . . . . . . . . . . 95 5.6 Position and Velocity proles for CHW Model (2.572 orbital period) . . . . . . . . 95 6.1 Principle 1 & 2 of NEF in action. The life-cycle of the input signal how its encoded and decoded back to originality . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 6.2 Approximation of dierent input signals . . . . . . . . . . . . . . . . . . . . . . . . 107 6.3 Approximation of dierent input signals . . . . . . . . . . . . . . . . . . . . . . . . 108 6.4 Block diagram for an LTI system. The integrator is driven by the signal _ x(t). Reproduced from Voelker and Eliasmith . . . . . . . . . . . . . . . . . . . . . . . . 109 6.5 Block diagram for an LTI system, with the integrator replaced by a rst-order lowpass lter. Reproduced from Voelker and Eliasmith . . . . . . . . . . . . . . . . 109 6.6 Adaptive PI controller for Inverted Pendulum system in Nengo GUI . . . . . . . . 111 x Abstract There has been an increasing interest in on-orbit autonomous servicing and repair of satellites as well as controlled active debris removal (ADR) in the space industry recently. One of the most challenging tasks for servicing/repair as well as for ADR is the rendezvous and docking with a non-cooperative tumbling resident space object (RSO). The work presented here deals with the optimal trajectory planning for proximity maneuvering and close approach for a nanosatellite with a tumbling RSO under high delity J 2 and quadratic drag perturbation model. The eect of J 2 and quadratic drag is studied on the relative motion and its importance as a requirement is established for implementing proximity operation algorithms for closed loop attitude control and relative navigation. Current work implements a minimum-control-eort problem for a 3-D rendezvous to a tumbling object which considers a three-degree-of-freedom model that consists of six kinematics states, position and velocity, for both the chaser/deputy and target/chief spacecraft. The control problem features additional path constraints with relative motion dynamics per- tinent to the proximity space operations to match the perching state vector of the nanosatellite over a feature of interest of the RSO. The path to close approach, with a terminal constraint of a small but nite positive relative speed at contact, is also discussed. The consequences of limited thrust and nite attitude maneuver time are taken into account and their eects on the closed loop translation control of the nanosatellite are analyzed.Moreover, the research also elaborates on the homotopic nature of the optimal control trajectories for variations in dierent mission design parameters including inclination, maximum thrust, and altitude. xi A secondary control problem is formulated with LQR cost, this facilitates a full state feedback controller with user dened state and control penalty. Optimal trajectories are studied under the dened penalty matrices. Finally, research shows a generic implementation of neuromorphic controller that is powered by population of neurons that has the ability to potentially do gain scheduling in realtime. xii Chapter 1 Introduction 1.1 Background There has been an increasing interest in satellite on-orbit autonomous servicing in the space in- dustry recently. JAXA recently completed a technology demonstration mission ETS7 [35]. NASA did an autonomous rendezvous mission through the Demonstration for Autonomous Rendezvous Technology (DART) in 2005 [23], where the mission was not completed due to more than expected propellant usage during rendezvous maneuvering. Air Force Research Laboratory's (AFRL) Ex- perimental Satellite Systems-10 and 11 (XSS-10 and XSS-11) have been developed in order to show the ability for a small sat to autonomously plan and rendezvous with a passive and cooper- ative Resident Space Object (RSO) in Low Earth Orbit (LEO) [6]. In addition, Defense Advance Research Projects Agency's (DARPA) Orbital Express (OE) Advance Technology Demonstration Program validated the technology and techniques for on-orbit refueling and conguration of two satellites [7][9]. The existence of this programs demonstrates that there is a need for a robust and eective autonomous close proximity control algorithms for multiple spacecraft. The use of micro-satellites to inspect, service, repair, deorbit and refuel larger spacecrafts is a long-term goal. In order to perform on-orbit service, the servicing spacecraft has to rst rendezvous and dock with the satellite to be serviced in orbit. The tumbling of an uncontrolled resident space object 1 (RSO) or the rotational and transnational motion of an non-cooperating RSO present challenges to the docking operations. This challenge is further elevated when the dynamics of the environment acting on the spacecrafts is highly non-linear. 1.2 Motivations The relative motion dynamics pertinent to proximity space operation only considers Hill's equa- tions, also known as Clohessy-Wiltshire relative equations of motion. The reason being Hill's equations are very simple to implement and have been successfully used to describe the relative motion of two satellite during rendezvous maneuvers. Additionally, this linear equations have a close form solution that relates the position of the spacecrafts with respect to their initial condi- tions. Because of this, they can be solved analytically and provide a solution that is fairly simple and easy to understand. These solutions also allow for an intuitive sense of the relative motion of satellites in clusters. Hill's relative equations are also therefore used in the design of control laws since the most eective control schemes require a set of constant coecient linearized dierential equations. These reasons and many more make them the logical rst choice in describing relative spacecraft motion. Although, the assumption made by Hill's equations that is cited repeatedly as the main source of error is that the Earth is perfectly spherical. Because the Earth is not perfectly spherical, but rather an oblate spheroid and the dominant term of the series expansion of this eect is the J 2 term, modeling errors are introduced by the non-central forces. In the derivation of Hill's equations as seen in the later chapter, the Earth is considered spherical, and the J 2 disturbance in not incorporated. Thus, they do not capture the J 2 disturbance eects. To reduce the eect of the modeling errors that are present in Hill's equations, much research has been done with varying success to incorporate the eects of the J 2 potential. Other papers do not use them at all, citing the fact that they do not capture the motion of the spacecrafts correctly. 2 However, while Hill's equations have proved very useful, they have several signicant limita- tions. Since they are linear, some error is introduced into the solution. However, linearization errors are not only limiting factor. Hill's equations are derived under the assumption that the disturbance forces acting on the satellites are also negligible. For rendezvous missions in the low earth orbit, signicant errors can be caused in the relative position due to drag. More over, while studying the motion of nanosatellites which are becoming increasing popular, there is very strong requirement for studying the eect of drag on the spacecrafts motion due to its size as even a small disturbance could cause it to sway away from it predetermined motion path. Moreover, for application concerned with proximity operations and rendezvous with non cooperating debris, it also becomes more crucial to implement drag disturbances as a part of the relative motion as a sense of reliability, to account for all perturbations so as to not create more debris. Much work has been done to implement drag perturbations for equations of motion which are linear in velocity for one body. To model a more realistic case, this thesis expands the work by implementing drag disturbance which is quadratic in velocity for the relative motion. It appears that there is a need for a set of linearized equations that are as easy and useful as Hill's equations, but at the same time capture the eect of the J 2 disturbance and drag force acting on the spacecrafts. In this thesis, a set of linearized, constant coecient dierential equations will be developed that capture the J 2 and quadratic drag perturbation. The importance of these perturbations will be studies by doing trajectory propagation. Finally, it will be shown for a control problem, the eect on the control-solution due to the formulated high delity dynamical constraints with respect to the simplied Hill's equations. These new equations are similar to Hill's equations in form. The radial and in-track motion is still coupled, but are separate from the cross-track direction. They are easily solvable, and can be implemented readily for rendezvous and proximity operation application for designing optimal control laws that take into account the perturbation both from J 2 and quadratic drag. 3 1.3 Literature Review Majority of research pertaining to applications that focuses on rendezvous and proximity operation with an uncontrolled RSO can be divided into two main segments. Firstly, the study of the relative motion dynamics to better model the environmental perturbations. Secondly, the study of dierent control schemes with varying complexity that comes about in a form of dierent path, event, and inequality constraints. Current proximity operation and rendezvous mission plans for spacecraft under the active debri removal have highlighted the need for improved dynamical formulations of the relative motion of spacecraft. The Clohessy-Wiltshire (C-W) equations are the rst relative motion equations for the rendezvous of spacecrafts. Linearizing the relative motion around a circular reference orbit, the C-W equations express the relative motion in terms of initial conditions that are dene in the Cartesian coordinates. A major drawback of the C-W formulation is the diculty in solving for the relative motion under any arbitrary perturbations. To overcome this problem, a substantial number of studies have examined the concepts for modeling motion under several disturbance forces. Tschauner and Hempel [36] removed the restriction of circular orbits and developed the linear equations of motion for any orbit around a spherical Earth. Also rst introduced by Lauden [19], time varying linear equations of motion utilizing an elliptical reference orbit. The resulting dierential equations, while no longer time-invariant can be solved analytically. Schaub and Alfriend [31], presents an analytical technique to derive a class of orbits that are termed J 2 invariant orbits. Satellites placed into J 2 invariant orbits are not immune to the J 2 disturbance force, but instead two satellites placed into these orbits will have the same drift, and thus the cluster will remain together. The J 2 invariant orbits are created by using two rst order conditions that determine the correct dierences in semi-major axis, eccentricity, and inclination. 4 In another paper by Schaub and Alfriend [30], an impulsive control algorithm is developed using mean orbital elements. The change from osculating orbital elements to mean orbital elements seems to have two purposes. The rst is that small periodic variations are not captured by the control system and thus propellant is not wasted counteracting these variations. Instead only the secular drift in any of the orbital elements is corrected for. The second purpose for switching to mean orbital elements makes the impulsive maneuvers being applied only once or twice during an orbital period. Sedwick, Miller and Kong [33] analyzed an assortment of perturbation eects of a satellite cluster. These included a non-spherical Earth, atmospheric drag, solar radiation pressure, and magnetic eld interactions. This is done by using a non-dimensional approach to determine the scaling of each perturbation eect. Sedwick [33][34] also describes the linear combination of Hill's equation with J 2 eect. The equations derived by Schweighart and Sedwick [33][34] provides insight into the relative motion of satellites under the in uence of the J 2 potential. The period mismatch between the in-plane and crosstrack motion produces an eect called tumbling in which the spacecraft formation appears to tumble around the orbital angular momentum vector. The formulated J 2 model with corrections described in next chapter, draws heavily from the studies conducted by Schweighart and Sedwick. The models presented by them are expanded to study proximity operations and are integrated with drag eects to better represent the environmental perturbations. Humi and Carter [14] [15] [17] [16] presents a formulation to implement drag eects in relative motion that is linear in velocity. It also expands on the scope by implementing liner drag for non- circular or low-eccentricity elliptic reference orbits. It later presents eects of quadratic drag on orbits that are not highly eccentric but decays as a result of drag. Much inspiration is taken from the research here in order to formulate relative motion models that implement drag perturbations that are quadratic in velocity, as discussed in next chapter. 5 With all the research pertaining to relative motion models for dierent applications, they all seem to fall into three general categories: papers that utilize Hill's equations, papers that show Hill's equations do not capture other perturbation eects, and papers that use a non-linear method of capturing the perturbation eects. Moreover, methods that are developed as a part for implementing J 2 perturbation, relies heavily on representing problem in dierent coordinate frames to simplify the derivation thus inducing even more linearization errors. Additionally, not a lot of research has been done that puts both the J 2 perturbation and drag that is quadratic in velocity together to form a higher delity model that captures both the perturbations. By surveying the current research, there appears to be a need for a set of linearized equations that are easy to use, but at the same time are able to capture both the J 2 disturbance along with quadratic drag experienced by the spacecrafts. From a theoretical standpoint, the thesis elaborates on the previous work by Boyarko et al. [11][12][13], who studied the minimum-time and minimum-control problem for a rendezvous of a chaser satellite to a tumbling object that considers a full six-degree-of-freedom model of both chaser and target. The current work expands the scope by adding the eect of J 2 geopoten- tial disturbance and the quadratic drag perturbation to investigates its eects on the proximity operations. 1.4 Organization This thesis is organized as follows: Chapter [1] contains a summary of past and current research on this topic and an outline of several concepts employed in the thesis. It also goes over motivations that characterizes the research outline. Chapter [2] starts by describing all the coordinate systems used for formulating the dierent relative motion models to better understand the problem of maneuvering one spacecraft relative 6 to another. The rst model described is the Hill's models or better known as Clohessy-Wiltshire relative equations of motion, assumptions pertaining to the models are explained. The later sections, question these assumptions and elaborates further by implementing time-averaged J 2 with corrections models. Finally, it is discussed for the scope of the application this thesis focuses on, that it is also benecial to implement drag disturbance as a part of the relative motion. All the sections here uses the linear orbit theory in order to derive the implementation of dierent perturbations on the system of spacecrafts. The last section under this chapter provides a new high delity linear model that encompasses both time-average J 2 with corrections and quadratic drag model. This model is the baseline model over which control problems in later chapters are implemented. Chapter [3] focuses on preliminary studies that shows the importance of formulating the above models. The purpose of this is to facilitate a better dynamical constraint for dierent control schemes discussed later in the thesis. Relative motion trajectories are simulated in time for dier- ent models derived in the chapter above. In case of J 2 disturbances, also shown are the in uence of the orbital period and nodal drift corrections as stated in the above chapter. Trajectory simula- tions for the full linear model is also presented. Finally, all the models are validated by simulating the trajectories at 0 inclination. The rational for choosing this inclination for validation purposes in also discussed in detail. Chapter [4] formulates the control problem and discusses dierent constraints that are im- plemented as part of the minimum-control-solution. Baseline parameters for the optimization problem including the selection of potential target for the proximity operation problem are also discussed in detail. Minimum propellant solution for the full linear model is presented and the eect of the control solution on varying orbital parameters are discussed. The physics of the formulated model in the above chapter is once again validated by comparing the control proles and the solutions from the full linear model with disturbance to the Hill's models. 7 Chapter 2 Relative Motion Models The purpose of this chapter is to describe the dierent relative motion models for the use of relative motion analysis so as to better understand the problem of maneuvering one spacecraft relative to another, when they are in close proximity. The rst section denes the coordinate frames that are used throughout the thesis. The section following explores one of the most common relative motion model, Clohessy-Wiltshire equations or better known as Hill's rendezvous equations while laying out the underlying assumptions for the CW model. The later two sections questions these assumptions and elaborates further on the derivation for implementing J 2 and Drag disturbance forces. Finally a full linear model is presented that captures both the described perturbations which is considered as the baseline model for studying the close-proximity control problem, described in later chapter. 2.1 Coordinate Frames This section denes the four coordinates frames used for the development of the relative motion models. The rst coordinate system is an inertial frame (ECI) in which absolute position of the reference orbit (chief/target) is dened. The other three frames are dierent body-xed coordinate 8 systems which are used to describe the relative position of the spacecraft (deputy/chaser) with respect to the reference orbit (chief/target). 2.1.1 ECI - Earth Centered Inertial Coordinate Frame The Earth centered inertial frame is described by ^ X ^ Y ^ Z. The ^ X ^ Y plane coincides with the equatorial plane, and the ^ Z axis points at North Pole completing the right hand rule. Figure [2.1] illustrates the ECI frame along with a reference orbit. ˆX ˆY ˆZ ˆr ˆ i ˆθ Figure 2.1: ECI Coordinate Frame and ^ r ^ ^ i Spherical Body-Fixed Coordinate Frame 2.1.2 ^ r ^ ^ i - Spherical Coordinate Frame Figure [2.1] also describes the ^ r ^ ^ i body-xed spherical coordinate frame. The ^ r points in the radial direction, away from the center of earth. The ^ i is the azimuthal angle measured along the line of nodes. And nally ^ is the co-latitude measured from the ascending node which acts 9 as the pole of the spherical coordinate frame. This coordinate system is used to describe the gravitational acceleration and the J 2 potential due to the spherical earth. 2.1.3 LVLH - Local Vertical Local Horizontal Coordinate Frame Figure [2.2] illustrates the local vertical local horizontal body-xed coordinate frame. This coor- dinate system is the foundation to all the relative models described below. A proximity operation maneuver usually involves a target vehicle (described by a reference orbit), which is passive and non-maneuverable, along with a chaser vehicle, which is active and performs the maneuvers re- quired ti bring itself alongside the target [24]. ˆX ˆY ˆZ r ref ˆN ˆT Reference Orbit Figure 2.2: LVLH Body-Fixed Coordinate Frame This moving coordinate system has its origin at the target as show in gure [2.2], the ^ Raxis is directed along the outward radial direction, the ^ Naxis is normal to the orbital plane of the target spacecraft and lies in the direction of the orbital angular momentum vector. Finally ^ Taxis completes the right hand rule such that it is perpendicular to both ^ R and ^ N and points in the direction of the target spacecrafts instantaneous velocity vector, tangential to the target spacecrafts orbit. 10 2.1.4 ^ x ^ y ^ z - Curvilinear Coordinate Frame The nal coordinate frame is a body xed coordinate system that is very similar to the LVLH coordinate system. Figure [2.3] illustrates the curvilinear ^ x (radial) ^ y (along-track) ^ z (cross- track) axis system. The ^ x vector points in the radial direction, the ^ z vector is perpendicular to the orbital plane and points in the direction of the angular momentum vector. Finally, the ^ y vector completes movement. In this ^ x ^ y ^ z coordinate system, the specication is made that it is a curvilinear coordinate system. The ^ x vector remains unchanged, however the ^ y and the^ vector 'curves' around the orbit. ˆx ˆy ˆz Figure 2.3: ^ x ^ y ^ z Curvilinear Body-Fixed Coordinate Frame Unlike the Local Vertical, Local Horizontal (LVLH) body xed frame, implemented by the Hill's equations, the only dierence comes from the fact that the LVLH frame is not dened as a curvilinear system but as a rectangular. Conversion between both these body xed coordinate systems is fairly straightforward. The radial vector in each coordinate system completely coincides. The ^ y vector, ^ T and ^ coincide; and the ^ z vector, ^ N and the ^ i coincide. The only dierence comes from the fact that the LVLH frame is not dened as a curvilinear system but as a rectangular. 11 2.2 Clohessy-Wiltshire Relative Motion Model Following is the derivation of the CW equations of motion using the linear orbit theory [24]. Figure [2.4] describes the relative position x between two spacecrafts in the LVLH frame as described in Section 2.1.3. ˆX ˆY ˆZ r ref ˆN ˆT Reference Orbit r x Figure 2.4: x - relative position between two spacecrafts dened in LVLH Assuming a 2-body motion, the equation of motion for the target body (reference orbit) can be written as, r ref = r ref r 3 ref (2.1) Similarly, the equation of motion for the chaser (secondary) spacecraft including another forces can be written as, r = r r 3 + u (2.2) Assumingjr r ref j and u are very small relative to dominant eects, jr r ref j<< r ref and u<< r ref r 3 ref (2.3) 12 This assumption of small dominant eects will be expanded upon in the next section by adding the eect of J 2 and drag disturbances. After accounting the above assumptions in Equation [2.1] and [2.2], the relative position acceleration between two spacecrafts can be written as, x = r r ref = r 3 ref " r ref r 3 ref r r 3 # + u (2.4) In order to get the above equation of relative motion only in terms of the reference orbit position, an approximation of r r 3 can be derived as follows, r = r ref + x (2.5) ) r r 3 = r ref + x (r 2 ref + 2 r ref x + x 2 ) 3=2 = r ref + x r 3 ref 1 + 2r ref x r 2 ref + x 2 r 2 ref 3=2 (2.6) ) r r 3 = r ref + x r 3 ref 1 + 2 r ref x r 2 ref + x 2 r 2 ref ! 3=2 (2.7) Rewriting the above equation in form of , 1 + 2 r ref x r 2 ref + x 2 r 2 ref ! 3=2 = (1 +) 3=2 (2.8) Because is small, Taylor series expansion is used about = 0, (1 +) 3=2 = 1 3 2 + 15 4 2 +O N (2.9) Substituting back for the value of in the above Equation [2.9] gives, 1 + 2 r ref x r 2 ref + x 2 r 2 ref ! 3=2 = 1 3 2 1 + 2 r ref x r 2 ref + x 2 r 2 ref ! + 15 4 1 + 2 r ref x r 2 ref + x 2 r 2 ref ! 2 +O N (2.10) 13 The above equation can be simplied further to condense other higher order terms inO N 1 + 2 r ref x r 2 ref + x 2 r 2 ref ! 3=2 = 1 3 r ref x r 2 ref +O N (2.11) Substituting the above Equation [2.11] in [2.7] gives the nal approximation of r r 3 r r 3 = r ref + x r 3 ref 1 3 r ref x r 2 ref ! +O N (2.12) Finally putting the above approximation in Equation [2.4] gives relative motion only in terms of r ref , x = r 3 ref " r ref r 3 ref r ref + x r 3 ref 1 3 r ref x r 2 ref ! +O N # + u (2.13) Expanding the above equation and compressing the higher order terms in x 2 , the nal approx- imation of the relative motion between two spacecrafts in the inertial frame can be calculated as, x = r 3 ref x + 3 r ref x r ref r ref r ref + u (2.14) Using the vector derivative rule, we can convert the above relative acceleration in the LVLH body-xed frame ( ^ R ^ T ^ N), x I = x R + (2! R _ x R ) | {z } Coriolis Term + ( _ ! R x R ) | {z } Euler Term +! R (! R x) | {z } Centrifugal Term (2.15) For a circular orbit assumption for the reference orbit. The angular velocity and the angular acceleration of the target spacecraft is, ! R =n ^ N and _ ! R = 0 (2.16) 14 wheren is the mean motion of the reference orbit. The Coriolis acceleration term and the centrifu- gal acceleration terms can be calculated as shown in Equation [2.20] and the Euler acceleration will go to zero given the uniform circular orbit assumption. The relative position x in the LVLH can be written as x R =x ^ R +y ^ T +z ^ N (2.17) ) _ x R = _ x ^ R + _ y ^ T + _ z ^ N (2.18) ) x R = x ^ R + y ^ T + z ^ N (2.19) With the formulation above, the Coriolis and centrifugal terms in the Equation [2.15] can be calculates as, 2! R _ x R = ^ R ^ T ^ N 0 0 2n _ x _ y _ z = 2n _ x ^ T 2n _ y ^ R (2.20) ! R (! R x) =! R ^ R ^ T ^ N 0 0 n x y z = ^ R ^ T ^ N 0 0 n ny nx 0 =n 2 x ^ Rn 2 y ^ T (2.21) Finally, putting Equation [2.20], [2.21], [2.14] and [2.15] together, the relative motion can be written as, r 3 ref 2 6 6 6 4 x + 3 r ref x r ref | {z } 3x r ref r ref |{z} ^ R 3 7 7 7 5 + u = x R + (2n _ x ^ T 2n _ y ^ R) + (n 2 x ^ Rn 2 y ^ T) (2.22) 15 Bringing x R on the left side and all the terms on the other side, and simplifying the relative motion term into its vector constituents. The above equation can be written as, x 2n _ y 3n 2 x =u x (2.23) y + 2n _ x =u y (2.24) z +n 2 z =u z (2.25) Here u = [u x u y u z ] T is the external force described in the LVLH frame. For a force free Hill's equation of motion u = 0 31 , and these constituents can be expressed as 3 scalar dierential equations, x = 2n _ y + 3n 2 x (2.26) y =2n _ x (2.27) z =n 2 z (2.28) The above model is used as the foundation for simulating the dynamics of the relative motion between two spacecraft for a verity of applications. Two fundamental aws that encompass this model is its underlying assumption of it only being valid for circular orbits. Secondly, the diculty of solving the equations under arbitrary perturbations. Moreover, the secular terms that come about as a part of the force-free solution cause non-periodic eects and will cause a violation of small relative orbit approximations. It is, therefore, an interest to evaluate alternative models that capture other perturbations better and can be applied to a more general set of orbits, so as to better model the dynamics for dierent applications. 16 2.3 High Fidelity J 2 Relative Motion Model The derivation begins with the analytic equation of motion for an orbiting satellite under the in uence of the J 2 potential. We use a similar coordinate system used to derive Hill's relative motion, the only dierence is replacing the ^ R ^ T ^ N rectangular frame with a ^ x^ y^ z curvilinear as shown is Figure [2.3]. r = g(r) + J 2 (r) (2.29) where, g(r) is the gravitational acceleration due to spherical Earth, g(r) =(=r 2 )^ r (2.30) J 2 (r) is the acceleration due to J 2 potential [24], J 2 (r) =(3=2)(J 2 R 2 e =r 4 ) 2 6 6 6 6 6 6 4 (1 3 sin 2 i sin 2 ) (2 sin 2 i sin cos) (2 sini cosi sin) 3 7 7 7 7 7 7 5 T [^ x ^ y ^ z] T (2.31) r is the position vector of the spacecraft, and ^ x-^ y-^ z is the coordinate system described in Figure [2.3]. The models are developed by considering the relative motion between one spacecraft and a reference orbit. For now, the derivation only requires that the reference orbit is constant radius (circular) and the equation of motion for the reference orbit is given by the Newton's law of motion in a central force eld: r ref = g(r ref ) (2.32) 17 Furthermore, linearizing the gravitational terms in Equation [2.1] with respect to the dened reference orbit results in r = g(r ref ) + g(r ref ) x + J 2 (r ref ) +rJ 2 (r ref ) x (2.33) where, the relative position of the satellite with respect to the reference orbit is given by x as represented in Figure [2.4] x = r r ref (2.34) When a spherical coordinate system (^ r - ^ - ^ i) with the pole aligned with the ascending node is used, the gradient of the g(r) gravitational acceleration can be calculated. The result is: rg(r) = 2 6 6 6 6 6 6 4 2(=r 3 ) 0 0 0 (=r 3 ) 0 0 0 (=r 3 ) 3 7 7 7 7 7 7 5 (2.35) The J 2 disturbance in Equation [2.31] is given in ^ x - ^ y - ^ z coordinates. However, the equation can be transformed directly to the ^ r - ^ - ^ i coordinate system without any loss of generality. Taking the gradient gives Equation [2.36] rJ 2 (r) = 6J 2 R 2 e r 5 2 6 6 6 6 6 6 4 (1 3 sin 2 i sin 2 ) sin 2 i sin 2 sin 2i sin 2 sin 2 i sin 2 1 2 sin 2 i( 1 2 7 4 sin 2 ) sin 2i cos 4 sin 2i sin sin 2i cos 4 3 4 + sin 2 i( 1 2 + 5 4 sin 2 ) 3 7 7 7 7 7 7 5 (2.36) Because the reference frame is rotating, rotational terms are added when calculating the relative acceleration and velocities of the satellite w.r.t. the reference orbit. x R = r r ref | {z } x I 2! R _ x R _ ! R x R ! R (! R x R ) (2.37) 18 For a circular reference orbit the rotation rate of the coordinate system is constant and given below. The Euler acceleration term can be ignored as a result of this assumption. ! R =n^ z; n = q =r 3 ref (2.38) Substituting Equation [2.33] into Equation [2.37] and rearranging the terms gives: x R + 2! R _ x R +! R (! R x R ) = g(r ref ) | {z } r ref +rg(r ref ) x + J 2 (r ref ) +rJ 2 (r ref ) x r ref (2.39) The equation above can be further simplied by canceling out g(r ref ) with r ref as they are equal as seen in Equation [2.32]. The following subsections now expand on the drawbacks of the above model. One major issue is the equations dependence on the gradients that are time varying. This drawback and other corrections are studied below and approximate solutions are derived so as to model the dynamics under J 2 perturbations even better with higher delity. Additionally, the following chapter 3 will show a comparison on the eect of modeling J 2 with and without these corrections on the relative position and velocity of the spacecrafts, to truly understand the importance of these corrections not just on the optimal solution but also for studying spacecraft trajectory propagation. 2.3.1 Time AveragedrJ 2 Equation [2.39] is a linearized dierential equation of motion with time varying coecients. The problem comes about thatrJ 2 term is not going to be constant except for equatorial orbits, where the eect of J 2 perturbation is negligible. A solution to this problem is to approximate the gradient of J 2 with the average value of the gradient term as a function of inclination for an entire period. Below is the approximate value of time averagerJ 2 term. 19 1 2 Z 2 0 rJ 2 (r)d = r 3 2 6 6 6 6 6 6 4 4s 0 0 0 s 0 0 0 3s 3 7 7 7 7 7 7 5 (2.40) here s is given by, s = 3J 2 R 2 e 8r 2 (1 + 3 cos 2i) (2.41) Adding the above time averaged J 2 term and removing g(r ref ) with r ref . The updated Equation [2.39] now becomes, x R + 2! R _ x R +! R (! R x R ) =rg(r ref ) x + J 2 (r ref ) + 1 2 Z 2 0 rJ 2 (r ref ) xd (2.42) By adding the time average of the gradient of the J 2 disturbance, the periodic component of the gradient is lost. As the values of g(r ref ), J 2 (r ref ) and the gradients of these terms are evaluated at the reference orbit, x must remain relatively small for the the rst order term of the Taylor series expansion and therefore, the mean radii and orbital rates of the perturbed spacecraft and the reference satellite must be similar. The fact that therJ 2 term has been time averaged over a period in Equation [2.40] has the eect to cause errors in cross-track motion. It is seen that cross-track motion is due solely to that the spacecraft orbit and the associated reference satellite are not co-planer. It is a periodic motion that is equal to zero when the two orbital planes intersect and is at a maximum 90 away from the intersection of the planes. As this research primarily deals with proximity and docking applications, where the satellite is in close vicinity with the reference orbit, both having similar inclinations, eects to the cross-track errors are not of interest as they will not be as signicant. Although addition and studying the eect of these cross track errors can be a topic that could be explored in future. 20 2.3.2 Correcting the orbital period of the Reference orbit Under the in uence of the J 2 disturbance, the perturbed spacecraft will have a dierent orbital period than when unperturbed. Because of this discrepancy, the perturbed spacecraft drifts from the unperturbed reference orbit and eventually the linearized equations breaks down. To x this problem, the period of the reference orbit must be adjusted to match the period of the spacecraft that is perturbed by J 2 . The change in period due to the J 2 disturbance as given by Schweighart and Sedwick [34] can be found from the average J 2 disturbance (not to be confused with the time average of the gradient of the J 2 term discussed in the previous subsection). The equation of motion of the reference orbit, as seen in Equation [2.32], now becomes, r ref = g(r ref ) + 1 2 Z 2 0 J 2 (r ref )d (2.43) where the time averaged J 2 acceleration is given by Equation [2.44] below, 1 2 Z 2 0 J 2 (r)d = 2 6 6 6 6 6 6 4 n 2 rs 0 0 3 7 7 7 7 7 7 5 (2.44) Heres is the same as in Equation [2.40] and is given by Equation [2.41]. For completeness of this subsection, it is as follows, s = 3J 2 R 2 e 8r 2 (1 + 3 cos 2i) (2.45) This results of averaging J 2 for a period is a non-zero value for acceleration in the radial ^ x direction only. This can be visualized as an additional force acting to increase the Kepelerian gravity term. Therefore, if a satellite is to remain in a circular orbit its orbital rate must be increased above that for Keplerian dynamics. 21 Substituting Equation [2.43] as another factor on the left hand side in Equation [2.42], x R + 2! R _ x R +! R (! R x R ) =rg(r ref ) x + J 2 (r ref ) + 1 2 Z 2 0 rJ 2 (r ref ) xd 1 2 Z 2 0 J 2 (r ref )d (2.46) As the rate of the reference orbit has changed, so has the average angular speed of the reference satellite, and the coordinate system which is based there. The new angular velocity can be found by equating the accelerations to give ! ! ! =nc^ z; c = p 1 +s (2.47) where s is given by Equation [2.45]. 2.3.3 Correcting the Reference orbit for nodal drift Although the preceding equations of motion are a vast improvement over Hill's equations when incorporating the J 2 disturbance, more can be done. Even though the orbital period of the reference orbit has been adjusted to match the period of the satellite under the in uence of the J 2 potential, they still drift apart due to separation of the longitude of the ascending node. Schweighart and Sedwick worked out analytically a separation distance between the perturbed satellite and the reference orbit as r ref 4 sini over one period. Similarly, after two orbital periods, the satellite and the reference orbit will be separated by r ref 24 sini, and this process will continue to cause them to drift farther and farther apart. Although the equations of motion in their current form do capture the bulk of this motion, the increasing separation causes the linearization to break down after several periods. Because of this, the reference orbit must again be redesigned so that it tracks this secular motion. 22 This separation in longitude of ascending node, shows that only the normal component of the J 2 acceleration is responsible for the drift. Applying the normal component of the J 2 disturbance acceleration to the updated reference orbit (Equation [2.43]) gives, r ref = g(r ref ) + 1 2 Z 2 0 J 2 (r ref )d + [J 2 (r ref ) ^ z]^ z (2.48) Adding this new force on the reference orbit, will facilitate an opposite drift between the two spacecrafts in the longitude of ascending node, so as to cancel the actual nodal drift, nally making the spacecrafts not being able to drift apart through time. As we add a component of the J 2 potential to the reference orbit for the correction of nodal drift, the component must be subtracted from the relative motion of one spacecraft with respect to the reference orbit. Equation [2.46] is now updated to, x R + 2! R _ x R +! R (! R x R ) =rg(r ref ) x + J 2 (r ref ) + 1 2 Z 2 0 rJ 2 (r ref ) xd 1 2 Z 2 0 J 2 (r ref )d [J 2 (r ref ) ^ z]^ z (2.49) 2.3.4 Time Averaged J 2 motion model with corrections The above Equation [2.49] can now be fully expanded to get the dierential equations that describe the relative motion of two spacecrafts under the in uence of time averaged J 2 with corrections in both orbital period and nodal drifts. x 2(nc) _ y (5c 2 2)n 2 x =3n 2 J 2 (R 2 e =r ref ) [ 1 2 3 sin 2 i ref sin 2 (kt)=2 (1 + cos 2i ref )=8] y + 2(nc) _ x =3n 2 J 2 (R 2 e =2r ref ) sin 2 i ref sin(2kt) (2.50) z + (3c 2 2)n 2 z = 0 23 where k is given below with J 2 = 1:081875 10 3 [?] and c is given in Equation [2.47], k =nc + 3 p J 2 R 2 e 2r 7=2 ref cos 2 i (2.51) As already mentioned under Section [2.3.1] the average of the dierential in J 2 is also likely to cause some cross-track errors, but since the magnitude of this error is primarily governed by the magnitude of the dierence in the inclinations between two spacecrafts. It is therefore, assumed to be very small for proximity operation applications which will have the inclination of both the spacecraft and the reference orbit be very similar. The above dierential equations are used for solving the optimal control problem described in the later chapter. These equations are also coupled with the quadratic drag relative motion equations to formulate the nal linearized equations of motion shown in Section [2.5]. Prior to setting up the simulation, next chapter goes over the eect of J 2 , with and without corrections as compared to Hill's equations, to highlight the motivation of this study and understand the eect of the errors on the relative position of the perturbed satellite. The next section now formalizes the eect of Drag into relative motion. For rendezvous and proximity operation applications with a nano satellite as a secondary spacecraft, it becomes crucial to analyze the drag perturbation as it could become one of the dominant force perturbing the spacecraft when considering it relative to its size. 2.4 Quadratic Drag Relative Motion Model As elaborated under Section [1.3] there already exist many models [14][15] [17] that integrate drag that is linear in the velocity. This section expands on the derivation on integrating the drag that is quadratic in the magnitude of the orbital velocity, as the assumption of drag that is liner in velocity is not as realistic. 24 Under the inertial frame discussed in Section [2.1.1], the equation of motion of the reference spacecraft is given by, r ref = g(r ref ) + f(r ref ) (2.52) In the above Equation [2.51], g(r ref ) is the 2-body gravitational acceleration resulting from a central force eld, and f(r ref ) describes the drag acceleration. The drag acceleration here is given in Equation [2.52] below, f(r) =f(r)j_ rj_ r (2.53) Here, is a physical constant associated with the drag coecient and the geometry of the space- craft. The functionf(r) represents atmospheric density and is assumed to only depend on altitude, its dependence on other variables like longitude, latitude, solar cycle are not been considered. As r represents the position of the reference spacecraft in the inertial frame, _ r is the velocity of the spacecraft such thatj_ rj = (_ r _ r) 1=2 . Substituting Equation [2.52] in [2.51] and rewriting r ref in terms of the the reference spacecrafts position and velocity, gives, r ref =g(r ref ) r ref f(r ref )j_ r ref j_ r ref (2.54) The functionsg andf here are continuously dierentiable and represents Newtonian Gravitational Field and the inverse law modeling the atmospheric density respectively. The function g is given by, g(r) = r 3 (2.55) and the function f is given by, f(r) = 1 r (2.56) 25 We can further simplify Equation [2.53] by writing the equation of motion in polar coordinates. The radial component is in the direction of r ref and the polar angle will be . From the angular momentum equation, the magnitude of the velocity is given by [40], j_ r ref j = ( _ r 2 ref + _ r 2 ref _ 2 ) 1 2 =r ref _ 2 6 6 6 6 4 1 + _ r ref r ref _ ! 2 | {z } O 3 7 7 7 7 5 1 2 (2.57) With an assumption of a reference orbit initially not highly eccentric. For circular or near circular orbits the magnitude of the radial velocity will be much smaller compared to the magnitude of the transverse velocity, i.ej _ rj << r _ . Using this assumption, Equation [2.56] can be linearized such that higher order terms under ( _ r) n =O can be ignored. Substituting Equation [2.56] with this linearized assumption into Equation [2.53], gives r ref =g(r ref ) r ref f(r ref )r ref _ _ r ref (2.58) Similarly the motion of the secondary (chaser/deputy) spacecraft around the reference satellite is given by, r =g(r)r f(r)j_ rj_ r (2.59) Here r is the position of the secondary spacecraft relative to the central force eld in the iner- tial frame as shown in Figure [2.4] and similar to depends on the drag coecient and the geometrical parameters of the secondary spacecraft. Now, the relative position of the secondary spacecraft relative to the reference in the inertial frame x is given by, r = x + r ref (2.60) 26 Substituting the above relation of spacecraft positions in the Equation [2.58], gives x + r ref =g(jx + r ref j) (x + r ref ) f(jx + r ref j)j _ x + _ r ref j ( _ x + _ r ref ) (2.61) The following assumptions can be used to simply the above model so as to linearize the relative motion between two spacecraft under quadratic drag. 1) As the secondary spacecraft is considered to be in close vicinity of the reference satellite for rendezvous and proximity operation applicationjxj<<jr ref j, any terms higher than the rst-order forjxj=jr ref j can be neglected 2) As previously discussed for the assumption of a reference orbit initially not highly eccentric, j _ rj<<r _ such that higher than rst-order term in _ r=r _ can also be neglected Accounting the above assumptions into the Taylor's expansion forg(jx +r ref j) andf(jx +r ref j) at r ref is given by, g(jx + r ref j) =g([(x + r ref ) (x + r ref )] 1 2 ) =g(r ref ) + dg(r ref ) dr ref r ref x r ref + jxj 2 2r ref +O N (2.62) Similarly, expanding the terms for f(jx + r ref j) gives, f(jx + r ref j) =f(r ref ) + df(r ref ) dr ref r ref x r ref + jxj 2 2r ref +O N (2.63) Substituting the above expressions in Equation [2.61] and neglecting all the higher order terms in jxj=jr ref j as per the assumptions above, gives x I =g(r ref )x + g 0 (r ref )[(r ref x)=r ref ]r ref ()f(r ref )r ref _ _ r ref f(r ref )r ref _ _ x f(r ref )[(r ref x)=r ref ] _ _ r ref f 0 (r ref )(r ref x) _ _ r ref (2.64) 27 Here g 0 () and f 0 () represents dierentiation with respect to r ref . It is to note that the above relative motion model is described in an inertial frame. In a body-xed rotating coordinate system as described in Section [2.1.3], the relative motion becomes, x I = x R + 2! R _ x R +! R (! R x R ) + _ ! R x (2.65) Removing the Euler acceleration and substituting Equation [2.64] in [2.65] while also applying the eects on rotation on the RHS which comes as a result of the derivative rule, gives, x R + 2! R _ x R +! R (! R x R ) =g(r ref )x + g 0 (r ref )[(r ref x)=r ref ]r ref ()f(r ref )r ref _ (_ r ref +! R r ref ) f(r ref )r ref _ ( _ x +! R x) [f(r ref )=r ref +f 0 (r ref )](r ref x)(_ r ref +! R r ref ) _ (2.66) Above equations can be written in the LVLH reference frame as a set of dierential equations, the ^ x radial axis points in the direction of r ref , ^ y is tangential to the orbit and point in the direction of the velocity vector, nally ^ z completes the coordinate system aligning itself with the angular momentum vector, normal to the orbital plane. x + f(r ref )r ref ! _ x = ! 2 +g(r ref ) +g 0 (r ref )r ref f(r ref ) _ r ref !f 0 (r ref )r ref _ r ref ! x + f(r ref )r ref ! 2 y + 2! _ y ()f(r ref )r ref _ r ref ! y + f(r ref )r ref ! _ y = ! 2 +g(r ref ) y 2! _ x 2f(r ref )r ref ! 2 +f 0 (r ref )r 2 ref ! 2 x ()f(r ref )r 2 ref ! 2 z +f(r ref )r ref ! _ z =g(r ref )z (2.67) 28 Finally, substituting for g(r ref ) and f(r ref ) from Equations [2.55] & [2.66] respectively, along with its derivatives with respect to r ref , the above dierential expressions take form, x = 2! _ y + 3h 2 r ref e 2 ! 2 x () ! 2 y +! _ x + _ r ref r ref !x +h 1 2 r ref _ r ref e ! y =2! _ x + () ! 2 x! _ y _ r ref r ref !yh 1 2 r 2 ref e ! 2 (2.68) z =! 2 z () ! _ z + _ r ref r ref !z The above transformation from Equation [2.67] to [2.68] is made possible by the close form solution developed by Carter and Humi [16] that describes the motion of the reference spacecraft under quadratic drag. The solution gives an orbit degradation equation of a spacecraft which is initially elliptical but not highly eccentric. r ref in the above equation will be given by this solution, r ref = h 2 (1 + 4 2 )= e 2 + 0 cos( 0 ) (2.69) Dierentiating the above equation with respect to time, we get _ r ref = (h 2 (1 + 4 2 )=) (e 2 2! 0 ! sin( 0 )) (e 2 + 0 cos( 0 )) 2 (2.70) with ! =hr 2 ref e as the angular velocity of the rotating coordinate system which is corrected in the presence of drag. h is the angular momentum that comes about as the constant of integral when deriving the orbit equation given in Equation [2.69]. While, 0 and 0 represents the initial eccentricity and the initial true anomaly of the reference orbit respectively. This model now lls the void in the literature by studying the drag perturbations that vary with the square of the magnitude of the velocity. A better approximation to this model is improving on the density modeling function f(r), the current density variation is captured by the inverse law which can be improved upon by considering atmospheric density decay that is exponential 29 with an increase in r. Moreover, a parametric curve t on the standard atmospheric data could also be considered . The diculty with more realistic density models is to be able to nd any close form solutions. Appendix [A] goes over the validation for this model, by making the initial eccentricity zero and considering both the spacecrafts as point mass making = = 0, it can be shown that the above Equation [2.68] reduces down to the CW relative equations of motion discussed in Section [2.2]. The next section now takes the linearization concepts shown in the above sections to integrate J 2 and quadratic darg together. The model described below will be the driving model used for studying dierent controls problems described in later chapters. As a part of the control problem, the relative motion described below will also act at the primary dynamical constraints for the control scheme studied. The eect of dierent factor contributing to the solution of the control problem under these dynamical constraints is also studied in Chapter [5] where a more analytical formation is presented. Prior to that, the next chapter provides comparisons of this model with Hill's relative equations in regards to trajectory propagation, and the importance to study and implement these disturbances. 2.5 Full Linear Model with Corrected J 2 and Quadratic Drag Dynamics For a fully linearized relative equations of motion that capture the J 2 perturbation with all the corrections as described in Section [2.3], along with capturing the quadratic drag disturbances as discussed above, we start with formalizing the equations of motion of the secondary spacecraft under these perturbations. In the inertial discussed in Section [2.1.1], the equation of motion of the secondary (deputy/chaser) spacecraft is given by, 30 r = g(r) + J 2 (r) + f drag (r) (2.71) Here g(r) is the gravitational acceleration due to spherical earth as described in Equation [2.30], J 2 (r) is the acceleration due to J 2 potential which is given by Equation [2.31] and nally f drag (r) is the acceleration due to drag quadratic in velocity as formulated in Equation [2.53]. For the completeness of this section, they are described as follow, g(r) =(=r 2 )^ r (2.72) J 2 (r) =(3=2)(J 2 R 2 e =r 4 ) 2 6 6 6 6 6 6 4 (1 3 sin 2 i sin 2 ) (2 sin 2 i sin cos) (2 sini cosi sin) 3 7 7 7 7 7 7 5 T [^ x ^ y ^ z] T (2.73) f drag (r) =f(r)j_ rj_ r = 1 r j_ rj_ r (2.74) The vector dierence between the secondary spacecraft and the target spacecraft (reference orbit) is shown as the relative position vector x, x = r r ref (2.75) Substituting r = r ref + x in Equation [2.71], r = r ref + x = g(r ref + x) + J 2 (r ref + x) + f drag (r ref + x) (2.76) In order to not lose the eect of drag constant parameter for the secondary spacecraft while undergoing the linearization through series expansion, the drag acceleration denition in the above equation is modied such that, 31 r = r ref + x = g(r ref + x) + J 2 (r ref + x) +f drag (r ref + x) (2.77) Here f drag (r) =f drag (r) from modifying Equation [2.74]. In order to linearize the above equation, the vectors g(), J 2 () and f drag () are expanded in a Taylor's series about the reference spacecraft r ref , which gives, g(r ref + x) = g(r ref ) +rg(r ref ) x +O N (2.78) J 2 (r ref + x) = J 2 (r ref ) +rJ 2 (r ref ) x +O N (2.79) Because the vector function f drag (r ref + x) which depends on the derivative of the linearized quantity _ r ref at which its expanded. Converting the vector function so it depends only in the position, gives f drag (r ref + x) =f(r ref ) _ r ref _ r ref | {z } f drag (r ref ) +f(r ref ) _ r ref _ x +f(r ref ) _ r ref " r ref x r 2 ref # _ r ref +f 0 (r ref ) r ref x r ref _ r ref +O r ref x r ref (2.80) Here f(r) =1=r the simplistic density modeling function as seen in Equation [2.74] with the negative sign carried inside. Ignoring all the second- and higher-order terms under the assumption that for a relative positionjxj=jr ref j<< 1, and substituting the above set of Equations [2.78] [2.79] and [2.80] into Equation [2.77], gives, x = g(r ref ) +rg(r ref ) x + J 2 (r ref ) +rJ 2 (r ref ) x +f drag (r ref ) +f drag (r ref ) " r ref x r 2 ref # +f(r ref ) _ r ref _ x +f 0 (r ref ) r ref x r ref _ r ref r ref (2.81) The inertial acceleration of the reference spacecraft can now be described as follow, 32 r ref = g(r ref ) + 1 2 Z 2 0 J 2 (r ref )d | {z } Period Correction + [J 2 (r ref ) ^ z]^ z | {z } Drift Correction +f drag (r ref ) (2.82) The correction terms comes about the derivation as formulated under Section [2.3] and the drag acceleration uses the modied denition with being the drag constant for the target spacecraft. Substituting Equation [2.82] in [2.81] and canceling out the g(r ref ) term, gives, x I =rg(r ref ) x + J 2 (r ref ) +rJ 2 (r ref ) x + ()f drag (r ref ) +f drag (r ref ) " r ref x r 2 ref # +f(r ref ) _ r ref _ x +f 0 (r ref ) r ref x r ref _ r ref 1 2 Z 2 0 J 2 (r ref )d [J 2 (r ref ) ^ z]^ z (2.83) The gradient of the vector functions in the above expression are expressed using Equations [2.35] and [2.36]. The relative motion above can now be converted into the rotating frame using the transport theorem. x I = x R + 2! R _ x R + _ ! R x R +! R (! R x R ) (2.84) Substituting Equation [2.83] into [2.84] and correcting for the derivatives on the RHS gives, x R + 2! R _ x R + _ ! R x R +! R (! R x R ) =rg(r ref ) x + J 2 (r ref ) +rJ 2 (r ref ) x +()f drag (r ref ) +f drag (r ref ) " r ref x r 2 ref # +f(r ref ) _ r ref ( _ x +! R x) +f 0 (r ref ) r ref x r ref (_ r ref +! R r ref ) 1 2 Z 2 0 J 2 (r ref )d [J 2 (r ref ) ^ z]^ z (2.85) The above now represents liner dierential equation that describes the motion of the secondary spacecraft with respect to the target spacecraft (reference orbit). With the substitutions for the gradient, the above equation can be written out in vector components are follow, 33 x = 2! _ y + (5c 2 2)n 2 x 3n 2 J 2 R 2 e r ref 1 2 3 sin 2 r ref sin 2 (kt) 2 1 + cos 2i ref 8 () ! 2 y +! _ x + _ r ref r ref !x +h 1 2 r ref _ r ref e ! y =2! _ x 3n 2 J 2 R 2 e 2r ref sin 2 i ref sin(2kt) + () ! 2 x! _ y _ r ref r ref !yh 1 2 r 2 ref e ! 2 z =(3c 2 2)n 2 z () ! _ z + _ r ref r ref !z (2.86) where k is given below with J 2 = 1:081875 10 3 [40], k =n p 1 +se | {z } ! + 3 p J 2 R 2 e 2r 7=2 ref cos 2 i (2.87) The angular velocity of the reference orbit which was assumed to be a circular initially is corrected to match the period due to the added J 2 perturbation, additionally there will also be a decay term associated with the angular velocity due to the drag perturbation as shown below, ! =n p 1 +se s = 3J 2 R 2 e 8r 2 (1 + 3 cos 2i) (2.88) The above equations now facilitates the study of the relative motion of two spacecrafts under a high delity perturbation model which includes both time-averaged J 2 with corrections and quadratic drag. This model is the foundation for the control problem implemented in the later chapter. Next chapter shows the importance of implementing these disturbances by propagating orbits in time for dierent criteria, magnitude of error in relative position and velocity between dierent models will show the in uence of these perturbations. Moreover, physics of the above formulated model is veried by simulating trajectories at 0 inclination compared to the Hill's model, similar behavior in expected in relative position and velocity. 34 Chapter 3 Eects Of Perturbation The purpose of this chapter is to understand the importance of the implemented perturbations that are modeled in the above chapter. This is done to facilitate a better dynamical constraint for dierent control schemes discusses later in this thesis. The rst section studies the eect of J 2 perturbation on the relative motion between two spacecrafts propagated through time. These linearized relative equations of motion with and without corrections for J 2 disturbances is also studies. The results are further compared to the Hill's model to verify the errors in relative position and velocity and give a sense of how crucial it is to implement J 2 perturbation along with its corrections. The second section similar to rst, presents simulated results for relative motion, now under quadratic drag perturbation. Furthermore, a comparison of these results with the Hill's model is presented to understand the magnitude of relative position and velocity errors. Finally, the last section presents a validation of the full linear model formulated above. The validation is carried out by simulating relative trajectory with reference orbit that is at 0 incli- nation and in considered initially in a circular orbit with e 0 = 0 35 3.1 Eects of the J 2 Disturbance This section highlights the eect of J 2 perturbations on the relative motion of the perturbed spacecraft with respect to the reference orbit. Moreover, errors related to reference orbital period and nodal drift are highlighted. Finally, the eect of disturbance of J 2 potential is studied in comparison to Hill's equation of motion. Numerical simulations are performed in MATLAB using xed-step ODE-5 integrator to plot rela- tive motion trajectories. The orbital parameters chosen for the reference orbit where the target is placed are: i ref = 70 , the apogee and perigee altitude is given by h a = 919-km andh p = 902-km respectively. The rationale for choosing the above orbital parameters is discussed in detail in next Chapter [4] where the optimization problem for this study is formulated. As the thesis pri- marily focuses on formulating the control problem for close proximity operation, an upper stage debri Agena-D is selected as the target and thus the orbital parameters of the debri are used for trajectory propagation for studying dierent relative motion models. The following Figure [3.1] demonstrates the eect of the J 2 geopotential disturbances at 70 inclination. The purpose to acquire the J 2 eect on the chaser from both the equatorial and polar regions, a logical choice of starting from [-1, 0, 0] km is employed as the initial condition for the relative position, thus in order to attempt to only simulate the J 2 eect on the pure in-track motion. The chaser is assumed to be at rest initially. All the simulations below for this section are propagated for the time span of two orbital period of the reference orbit which comes to 3.44 hrs. It is seen in Figure [3.1] the considerable relative position error created by the J 2 disturbance and the importance for its implementation while studying the optimal trajectory design for ren- dezvous and proximity operations. Additionally, it is noted that the J 2 model without corrections corresponding to Equation [2.39] adds sizable cross-track errors as a result of orbital period mis- match between the target and the chaser orbits and, the nodal drift caused by the separation of the longitude of the ascending node. 36 Figure 3.1: Relative motion trajectories with i ref = 70 for 2-orbital period The relative position errors between the commonly used Hill's model and the time-averaged J 2 with corrections model is shown in Figure [3.2]. As seen in the gure, the error in the in-track motion (^ x) tops out at 15 km and error propagation in time looks periodic. Similarly, for the radial motion (^ y) the error due to the J 2 seems to be increasing sharply with time. Finally, for the cross-track motion (^ z) even though the error is not signicant compared to in-track and radial motion, the prole shows an increasing deviation with time. With periodic errors in the in-tack motion and the increasing errors in both the radial and cross-track direction, a signicant eect of J 2 perturbation on the chaser is seen. 37 Figure 3.2: Relative position errors between Hill's Model and Time-averaged J 2 with corrections Model (i ref = 70 and 2-orbital period) Figure [3.3] presents similar relative position errors between the time-averaged J 2 model with corrections and, the J 2 equations without any corrections (Eq. [2.39]) to highlight the eect caused by the orbital period mismatch and the separation in the ascending node causing the nodal drift. As expected, the position error in all three directions is signicant. As a consequence, the imple- mentation of the corrections is crucial for studying the optimal trajectory design for proximity applications. Moreover when compared to Figure [3.2], it is seen that the implementation of the 38 corrections reduces error in the ^ z direction which causes the motion of the chaser to create heavy cross-track deviations from the nominal as seen in Figure [3.1]. Figure 3.3: Relative position errors between Time-averaged J 2 with corrections Model and J 2 without corrections [Eq. 2.39] (i ref = 70 and 2-orbital period) One way to validate the physics for the above formulated models is to analyze them at 0 inclination. As the eect of J 2 is not prominent on the equatorial plane, the time-averaged J 2 model with corrections model should behave similar to Hill's equations of motion. Figure [3.4] shows the relative motion trajectories for 0 under the in uence of the J 2 distur- bance force. Is is seen like expected that the relative motion of the chaser under the time-averaged 39 Figure 3.4: Relative motion trajectories with i ref = 0 for 2-orbital period J 2 with corrections model resembles with that of Hill's equations, as logically there would not be any eect of J 2 across the equator. Although a signicant error is seen in the relative position given by the J 2 equations without corrections. Figure [3.5] and Figure [3.6] illustrates the errors in the relative position of the spacecrafts between J 2 model with corrections with Hill's Equations and J 2 model without corrections re- spectively. As seen in Figure [3.5], the position errors between J 2 corrections model when compared with Hill's relative motion model are very small. The insignicant error noticed corresponds to the curvilinear coordinated system used for the derivation of the time-averaged J 2 equations. Validating the time-averaged J 2 corrections model once again at 0 inclination proves that the model used for the optimization problem formulation in the following section is a generalized case of Hill's Equation with J 2 perturbation incorporated in it. Figure [3.6] further shows the inability of the J 2 model without corrections to fully capture the orbital period mismatch and nodal drifts errors. As expected the position errors in all three directions are fairly signicant. 40 Figure 3.5: Relative position errors between Hill's Model and Time-averaged J 2 with corrections Model (i ref = 0 and 2-orbital period) From the study conducted in this section we see the importance of implementation of J 2 in the linear equations of relative motion for the study of the optimization operations. The next subsection continues by studying the eect of drag that is quadratic in velocity, the goal is again to characterize this disturbance as one of the dominant perturbation that can cause signicant errors in the relative position and velocity for studying optimal control problems for proximity operations. 41 Figure 3.6: Relative position errors between Time-averaged J 2 with corrections Model and J 2 without corrections [Eq. 2.39] (i ref = 0 and 2-orbital period) 3.2 Eects of Drag Disturbance This section highlights the eect of quadratic drag on the relative motion of the perturbed space- craft with respect to the reference orbit (target spacecraft). Parameters for the numerical simula- tions are similar to the Section [3.1] described above. The reference orbit will have an inclination of 70 with the perigee and apogee altitude given by h p = 902-km and h p = 919-km respectively. 42 Additionally, while implementing drag eect it is also crucial to dene the geometric param- eters for both the spacecraft, we start with assuming the secondary spacecraft (chaser) with a mass of m chaser = 24-kg and surface area of A chaser = 0:23 0:23 = 0:0526-m 2 . The rational about this assumption is discussed in detail under Chapter [4] where the control problem is for- mulated. As discussed previously, the reference space object is assumed to a upper stage debri Agena-D, the dry mass for this upper stage is given as m ref = 673-kgs with the cross section area ofA ref = 1:5 6:5 = 9:75-m 2 [37]. The ballistic coecient can now be calculated and given below in Equation [3.1]. chaser = 0:0526 10 6 24 = 2:204 10 8 km 2 =kg ref = 9:75 10 6 673 = 1:449 10 7 km 2 =kg (3.1) From the ballistic coecients it can be noted that the reference space object is almost 6:57 bigger than the secondary spacecraft. Section below further expands by studying the in uence of drag for a system of two spacecraft with higher ballistic coecient ratios. To simulate the most extreme scenario, the drag coecient of C d = 2:1 is used, this implies both the spacecrafts being at plates and are aligned perpendicular exposing their full area in the direction of motion. The secondary spacecraft is assumed to be at rest initially and starting at [1; 0; 0] km. All simulations below for this section are propagated for the time span of two orbital period of the reference orbit which comes to 3:44 hrs. The following Figure [3.7] demonstrates the eect of Quadratic drag at an altitude of 911-km. As seen the relative position and velocity errors between the drag model and Hill's equation of motion are very fairly noticeable. It is expected that the magnitude in error when compared to J 2 model with corrections from the above Section [3.1] will be less signicant, as the eect of J 2 disturbance relative to drag is highly dominant. 43 Figure 3.7: Relative motion trajectories at 900-km altitude with i ref = 70 for 2-orbital period The relative position errors between the Hill's model and the quadratic drag model is shown in Figure [3.8]. It can be seen that the position errors in the cross-track (^ z) are negligible for the initial conditions chosen. It is shown later with dierent initial conditions, the cross-track motion errors are in-fact signicant and also periodic with secular eects that grows with time. The radial (^ y) and in-track (^ x) motion shows signicant errors in relative position. Both seem to increase sharply with time. With relative position errors in all three motion, a signicant eect of drag perturbation on the chaser is seen. It is also to note that for the J 2 perturbation model, the relative position by Hill's model were always understated, whereas for the quadratic drag model, the relative positions by Hill's model are exaggerated. As drag force can be characterized as a friction force action opposite on the motion of the spacecrafts, the above exaggerated position errors by Hill's model seem to validate that intuition. As previously stated the magnitude of the errors when compared to J 2 disturbance with corrections are not too big. But the eect of drag perturbation on the spacecrafts vary on multiple factors. Hence, it facilitates a more through study of the quadratic drag model at dierent altitudes and for longer periods. Later in the section, relative position errors between dierent altitudes and dierent times periods are compared. 44 Figure 3.8: Relative position errors between Quadratic Drag Model and Hill's Model (i ref = 70 with 900-km Altitude for 2-orbital period) The following Figure [3.11] shows the relative position errors with change in initial conditions. The starting point for the chaser is now [1; 0;1]. This highlights the cross track (^ z) motion errors. As previously stated the errors are both signicant and periodic in nature. As seen in the Figures [3.9] and [3.10] above, the errors in the relative motion shows signicant drift contributed by the drag perturbation. These drifts also seem to increase with each orbital 45 Figure 3.9: Relative motion trajectories - i ref = 70 with I.C. = [1; 0;1] for 2-orbital period (900-km altitude) Figure 3.10: ^ x ^ z view for I.C. = [1; 0;1] for 2-orbital period (900-km altitude) - Drift view 46 period of the reference orbit. Figure [3.11] further shows the relative position errors between the quadratic drag model and the Hill's model. Figure 3.11: Relative position errors between Quadratic Drag Model and Hill's Model for I.C. = [1; 0;1] (i ref = 70 with 900-km Altitude for 2-orbital period) As stated earlier, because the drag perturbations depends on a lot of factors, it is crucial to study the eect of varying altitude on the relative errors. The eect of drag should get more dominant with decrease in altitude. Similarly, from the position errors seen above, it would be 47 also benecial to see the in uence of drag for longer orbital periods. The following Figures [3.12] and [3.15] shows the relative motion for dierent altitudes and times respectively. Figure 3.12: Relative motion trajectories with Quadratic Drag Model at multiple altitudes for I.C. = [1; 0;1] (i ref = 70 for 2-orbital period) It can be seen in the Figures [3.12] and [3.13] the relative position deviations cause due to the change in the altitude. To highlight the periodic error prorogation in the cross track motion, the initial conditions for the secondary spacecraft are chosen as [1; 0;1]. The deviations seen veries the intuition that at lower altitudes there is a higher eect of drag perturbations. This eect is clearly highlighted in Figure [3.13] which zooms at the nal state of the spacecrafts after 2 orbital-period. It is expected that at higher altitude the position compared to lower altitude would be exaggerated. This is purely because at higher altitude their would be simply less drag due to lower atmospheric density. Also to note that the position deviations between lower altitudes (between 500-km and 300-km) is much higher than the position deviations 48 between higher altitudes (between 900-km and 500-km) as shown in Figure [3.14]. This can be again regarded to the fact that at higher altitudes the eect of drag will be much lower compared to lower altitudes. Figure 3.13: Zoomed view of Relative motion trajectories with Quadratic Drag Model at multiple altitudes for I.C. = [1; 0;1] (i ref = 70 for 2-orbital period) The following Figures [3.15] and [3.16] highlights the eect of longer periods under the drag perturbation. 49 Figure 3.14: Relative position errors with Quadratic Drag Model for dierent altitudes for I.C. = [1; 0;1] (i ref = 70 for 2-orbital period) 50 Figure 3.15: Relative motion trajectories - i ref = 70 with I.C. = [1; 0;1] for 7-orbital period (900-km altitude) As discussed previously, the errors in relative position increase sharply with time in both radial (^ y) and in-track (^ x) motion, this is validated in the Figures [3.15] [3.16] and [3.17]. When compared to Figure [3.10] that shows the drift for 2 orbital-period, it can be seen that for longer periods here the drift due to the drag perturbation becomes signicantly prominent with time. Once again, the errors in the cross track (^ z) motion are simply periodic in nature and does not have any secular growth associated with it in time. Majority of the errors in the relative position comes about the radial motion as seen in Figure [3.17]. From the study conducted in this section we see the importance of implementing quadratic drag in the linear dierential equations of relative motion for the study of optimal control problems. Similar to how the time-averaged J 2 model was validated in the above section, we can imple- ment a similar validation study to cross examine the correctness of the physics in the quadratic 51 Figure 3.16: ^ x ^ z view for I.C. = [1; 0;1] for 7-orbital period (900-km altitude) - Drift view Figure 3.17: Relative position errors between Quadratic Drag Model and Hill's Model for I.C. = [1; 0;1] (i ref = 70 with 900-km Altitude for 7-orbital period) 52 drag model. Unlike J 2 the eect of drag does not depend on the inclination of the reference orbit, but more so on the geometrical parameters of both the spacecrafts and the altitude chosen for the reference orbit as seen above. Therefore, in order to validate the quadratic drag model, the drag constants for both the spacecrafts as dened in Equation [2.59] can be set to zero i.e. = = 0. In this case, the only eect that will be pertaining to the relative position errors will be due to the fact that the Hill's models assumes a circular reference orbit. The position errors for this case are highlight in Figure [3.18] and [3.19]. Figure 3.18: Relative motion trajectories for I.C. = [1; 0; 0] (i ref = 70 with 900-km Altitude for 2-orbital period) and = = 0 As previously stated, the relative position errors seen in the above Figure [3.19] are caused due to the reference orbit being eccentric. This can be validated by further setting the apogee and perigee altitude equal to each other i.e. h a = h p = 902-km. For a case like this, is can 53 Figure 3.19: Relative position errors between Hill's Model and Quadratic Drag Model for I.C. = [1; 0; 0] (i ref = 70 with 900-km Altitude for 2-orbital period) and = = 0 be said the quadratic drag model should behave exactly like Hill's equations of motion for small orbital-periods. As the models formulated in the previous chapter are linearized in nature, there is a possibility that for longer periods, there will be some linearization errors that will become 54 dominant in time. For the case with smaller orbital-periods Figure [3.20] and [3.21] shows the relative position errors. Figure 3.20: Relative position errors between Hill's Model and Quadratic Drag Model for I.C. = [1; 0; 0] (i ref = 70 for 2-orbital period) with = = 0 and h a =h p = 902-km 55 3.3 Validating the Full Linear Model This section implements the full linear model as described in Equations [2.86]. As previously discussed the importance of characterizing J 2 and quadratic drag disturbances, this goal with this section is to propagate trajectories by including both the perturbations into the relative motion. Simulation parameters are similar to previous two sections. The reference orbit will have an inclination of 70 with the perigee and apogee altitude given by h p = 902-km and h a = 919-km respectively. Moreover, the geometric parameters for the spacecraft as similar to the one described under Section [3.2]. The secondary spacecraft is assumed to be at rest initially with a starting point at [1; 0; 0] km. All simulations below for this section are propagated for the time span of two orbital period of the reference orbit. Figure 3.21: Relative motion trajectories for Full Linear Model with both J 2 and quadratic drag perturbation for 2 orbital-period (i ref = 70 and 900-km altitude) with IC = [1; 0; 0] As seen in Figure [3.21] the eect of J 2 and drag on the relative motion is very signicant. When comparing the relative motion errors described in Figure [3.22] to Figure [3.2], it can be 56 shown that they are highly comparable. This comes about because of the dominant eect of J 2 disturbance when compared to quadratic drag at the given altitude. For the case when the inclination of the reference orbit is 0 , drag perturbation on the spacecraft will be much more dominant than the J 2 disturbance. Figure 3.22: Relative position errors between Hill's Model and Full Linear Model with both J 2 and quadratic drag perturbation for 2 orbital-period (i ref = 70 and 900-km altitude) with IC = [1; 0; 0] 57 We can similar to previous sections validate the physics of the full linear model, this can be done by simulating a trajectory for 0 inclination and setting the drag constant for both the spacecraft to be = = 0. If the drag constants are not set to zero, the eect of drag for the given altitude will be accounted for, giving similar results as demonstrated in the above Section [3.2]. In such a case as previously states, the drag eect would be more dominant compared to the J 2 perturbation. For the case of validation, Figures [3.23] and [3.24] shows the relative position errors between Hill's equations and the full linear model. Figure 3.23: Relative motion trajectories for Full Linear Model with both J 2 and quadratic drag perturbation for 2 orbital-period (i ref = 0 and 900-km altitude) with IC = [1; 0; 0] and = = 0 As seen in Figures [3.23] and [3.24] the relative position errors are signicantly smaller when compared to the Figure [3.22] which simulates eect of both J 2 and drag perturbations. The small errors as seen can be attested to the non circular reference orbit for the full linear model. These errors are also very similar to the errors shown in Figure [3.4] which are the validation errors of Hill's equation with respect to the time-averaged J 2 with corrections. 58 Figure 3.24: Relative position errors between Hill's Model and Full Linear Model with both J 2 and quadratic drag perturbation for 2 orbital-period (i ref = 0 and 900-km altitude) with IC = [1; 0; 0] and = = 0 With the preliminary studies conducted in this chapter, it becomes very evident for imple- menting the eects of both J 2 and quadratic drag perturbations into the relative motion models. For the purpose to facilitate a better dynamical constraint model for the study of optimal control laws. Next chapter takes the full linear model and proposes a minimum-propellant-solution for studying proximity operations. 59 3.4 Comparisons with Nonlinear Model This sections shows and validates the eectiveness of the developed full linear model. Nonlinear dynamics are fully integrated to compare the relative position and velocity errors between the nonlinear motions versus the developed linear model dynamics. Both short t = 2T 0 and longer periods t = 5T 0 trajectories are compared to realize the break in linearity for longer time span orbits. Simulation parameters are similar to previous two sections. The reference orbit will have an inclination of 70 with the perigee and apogee altitude given by h p = 902-km and h a = 919-km respectively. Moreover, the geometric parameters for the spacecraft as similar to the one described under Section [3.2]. The secondary spacecraft is assumed to be at rest initially with a starting point at [1; 0; 0] km Figure 3.25: Relative motion trajectories for Full Linear Model with both J 2 and quadratic drag perturbation for 2 orbital-period (i ref = 0 and 900-km altitude) with IC = [1; 0; 0] compared to Nonlinear Dynamics Model 60 Figure 3.26: Relative position errors between Nonlinear Dynamics Model and Full Linear Model with both J 2 and quadratic drag perturbation for 2 orbital-period (i ref = 0 and 900-km altitude) with IC = [1; 0; 0] As seen in Figures [3.25] and [3.26] full linear model tracks the nonlinear dynamics very well for short periods, relative position errors in the cross-track motion looks to be increasing with each period. These cross track errors can be corrected for longer periods, given the context of the research primarily revolves around rendezvous and docking operations which should not last for longer orbits, it was chosen to not correct for the cross track errors. As ever increasing the cross-track error trend seems, the relative errors in the cross track motion are still signicantly smaller in magnitude at around 0:02-meters. Figure [3.27] and [3.28] shows similar results for when the orbits are inclined. 61 Figure 3.27: Relative motion trajectories for Full Linear Model with both J 2 and quadratic drag perturbation for 2 orbital-period (i ref = 70 and 900-km altitude) with IC = [1; 0; 0] compared to Nonlinear Dynamics Model As seem in the above Figure [3.27] full linear model captures most of the non-linearity in the dynamics. Relative position errors seem bounded for shorter orbits at higher inclinations. Cross-track errors seem unbounded and increasing but as discussed previously these errors are signicantly small in magnitude compared to its radial & in-track counterparts as seen in Figure [3.28] Figure [3.29] shows eect of long period on the linearity. It is clearly seen for longer periods up until 5-periods at lower inclinations Full Linear Model seems to capture majority of the non-linear dynamcis. In-plane errors y error are the major contributing factor in the relative position errors as seen in Figure [3.30]. These large errors in the in-plane motion are the eect of time averaging the J 2 pertubations for 1 orbital period. The nonliner eect of J 2 at higher inclination is very dominant and overtakes the time averaging assumption. 62 Figure 3.28: Relative position errors between Nonlinear Dynamics Model and Full Linear Model with both J 2 and quadratic drag perturbation for 2 orbital-period (i ref = 70 and 900-km altitude) with IC = [1; 0; 0] 63 Figure 3.29: Relative motion trajectories for Full Linear Model with both J 2 and quadratic drag perturbation for 5 orbital-period (i ref = 0 and 900-km altitude) with IC = [1; 0; 0] compared to Nonlinear Dynamics Model 64 Figure 3.30: Relative position errors between Nonlinear Dynamics Model and Full Linear Model with both J 2 and quadratic drag perturbation for 5 orbital-period (i ref = 0 and 900-km altitude) with IC = [1; 0; 0] Figures [3.31] and [3.32] shows the inability to fully capture the non-linear dynamics at higher inclinations for longer orbital periods. For smaller inclinations Full Linear Model is still a valid model as the the maximum relative postion errors are around 6-meters. 65 Figure 3.31: Relative motion trajectories for Full Linear Model with both J 2 and quadratic drag perturbation for 5 orbital-period (i ref = 70 and 900-km altitude) with IC = [1; 0; 0] compared to Nonlinear Dynamics Model This section highlights the correct usage of the full linear model. Next chapters build on this and dierent control formulations are studied. Given that the developed linear models captures majority of the non-linear dynamics in the motion at smaller orbital periods, all the control formulations here on out will impliment optimal simulation bound to be less than 2 orbital periods. 66 Figure 3.32: Relative position errors between Nonlinear Dynamics Model and Full Linear Model with both J 2 and quadratic drag perturbation for 5 orbital-period (i ref = 70 and 900-km altitude) with IC = [1; 0; 0] 67 Chapter 4 Control Problem Formulation This Chapter develops a model of target-chaser rendezvous. The LVLH (^ x ^ y ^ z) coordinate system described in Section [2.1.3] is used for the model formation. We start from the arbitrary relative position and would like to bring the two spacecraft together for proximity operation. Using the formulations described in Chapter [2] based on the full linear model which incorpo- rates both time-averaged J 2 with corrections and quadratic drag perturbations, the dynamics of the two system can now be described as follow. The transnational kinematics and dynamics of a chaser spacecraft in the orbit frame centered on the target vehicle (reference orbit) are given by Equation [2.86], x = 2! _ y + (5c 2 2)n 2 x 3n 2 J 2 R 2 e r ref 1 2 3 sin 2 r ref sin 2 (kt) 2 1 + cos 2i ref 8 () ! 2 y +! _ x + _ r ref r ref !x +h 1 2 r ref _ r ref e ! + f x y =2! _ x 3n 2 J 2 R 2 e 2r ref sin 2 i ref sin(2kt) + () ! 2 x! _ y _ r ref r ref !yh 1 2 r 2 ref e ! 2 + f y z =(3c 2 2)n 2 z () ! _ z + _ r ref r ref !z +f z (4.1) where f x ;f y and f z are applied forces (controls) expressed in the LVLH coordinate system. 68 Equation [4.1] hence denes a 6-state system of dierential equations of relative motion gov- erning the control problem dynamics. Combined into the state vector X these states are, X = [x;y;z; _ x; _ y; _ z] T (4.2) The governing dynamics assumes three normalized controls: u = f x f xmax ; f y f ymax ; f z f zmax (4.3) For simplicity, it is assumed that f imax = 1m=s 2 for i = x;y;z. Once again, these controls are three normalized components of a transnational force acting on a chaser f i (i =x;y;z), expressed in LVLH coordinate frame. All three controls are bounded: u min u u max . The details of the thrust bounds are discussed in the Subsection [4.0.1] along with chaser mass and its thrust parameters. Using the controls the two spacecrafts are brought together from some initial condition, given by 6 initial value of states x i (t 0 );i = 1;:::; 6 to end point constraint described by a nal state, which for this problem is assumed to be a null i.e. X f = 0. Figure 4.1: Space vehicle thruster conguration as l 1 norm As the spacecrafts control and guidance strategies are primarily dictated by the amount of pro- pellant the proposed control scheme is required to accomplish the specied mission, the condition of the objective function therefore becomes propellant minimization. The way this is achieved in 69 this thesis is by minimizing the total thrust acceleration under theL 1 norm of the control. Thrust- ing here is achieved by six (ungimbled) identical engines rigidly mounted to the body axes as seen in Figure [2.1]. The above is used as the thruster conguration for the secondary spacecraft. Equation [4.4] now describes the cost functional for the minimization of thrust acceleration. J(f; u) = Z t f t0 (ju 1 (t)j +ju 2 (t)j +ju 3 (t)j)dt (4.4) 4.0.1 Optimization Problem Setup A sample maneuvering scenario is considered with the chaser center of mass starting at a distance of 1-km behind, 1-km below, and 1-km sideways from the target center of mass. Hence the initial values of the state for computer simulations discussed in the following sections is: X(t 0 ) = [1000;1000;1000; 0; 0; 0] T (4.5) The endpoint constraint is state constraint at the nal time. For simplicity it is assumed that the nal state is at the CoM of the target and with null relative velocity vector. X(t f ) = [0; 0; 0; 0; 0; 0] T (4.6) Baseline parameters are selected from a proposed controlled active debris removal (ADR) mission Curimba by Udrea and Nayak [37]. The authors describes Agena-D RBs as one of the potential target for the mission. Hence, the orbital parameters for this upper stage debri corresponds to the reference orbit chosen for this problem. As most of the Agena-D RBs are populated at approximately 70 the reference inclination for this problem is assumed 70 . Additionally, the apogee and the perigee altitude selected for the reference orbit where our target spacecraft is placed resembles that of the Agena-D RB which corresponds to r a = 919-km and r p = 902-km 70 respectively. Furthermore, a time constraint is added to the simulation by xing the nal time to one-orbital period of the reference orbit, as for the ADR mission it requires the chaser to rendezvous with the target with the shortest amount of time. The mission proposes a preliminary design of the chaser by considering a 12U CubeSat. This CubeSat is a three-axis stabilized spacecraft of 240240360 mm and a mass of 24 kg. It is known that the CubeSat has a miniaturized system consisting of 16 solenoid-fed thrusters placed such that they generate positive and negative torques about all three body axes with dual redundancy. The thruster arrangement also generates positive and negative forces about two of the three body axes. For the purpose of this work it has been assumed that the force along the third body axis is available when required. Future research will include constraints on the generation of force about the third body axis to take into account the nite slew time. It is known from the preliminary design provided by the authors that the propellant used will be 1,1- di uorethane which can generate about 70 s of specic impulse. As for the rst mission design iteration they assumed that each thruster produces 50 mN of thrust force [37]. Knowing now the thrust force for individual thrusters (50 mN) and the chaser mass (24 kg), control bounds can be calculated for the optimal problem for the Equation [4.3]. By these means the control bounds thus will be calculated byF T =m chaser , where F T is the total thrust produced by the chaser spacecraft thrusters. Solving for the control bounds for the dened rendezvous scenario and the assumed chaser parameters gives us:8:33 mm/s 2 u 8:33 mm/s 2 The above described optimal control problem is solved numerically, by the Gauss Pseudospec- tral Optimization Solver (GPOPS) [25]. GPOPS is an open source code that implements a direct collocation method based on the Gauss Pseudospectral approach for solving optimal control prob- lems. It also employed hp-adaptive methods to better approximate the solution by increasing both the order of the tted polynomial and by increasing the collocation points when needed. Propellant-optimal trajectories are simulated with full linear model described in Equation [4.1]. To validate the correctness of the derived model, trajectories with 0 inclination are simulated 71 and compared with Hill's equations along with its control proles. It should be seen that for this specic test condition, the optimal trajectories in both case should be perfectly similar along with its trust prole. Moreover, the eects of the correction errors from the J 2 perturbation model on the optimal trajectories is also shown. The simulation results obtained including the optimal control proles are discussed below. 4.1 Minimum Control Solution for Full Linear Model The minimum-propellant-control solution with the optimal controls dened in Equation [4.1] is here presented rst. Additionally, the eect of the optimal trajectories with respect to the inclination is discussed, it also shows the eect of the full linear model at dierent inclinations and how the optimal trajectories are in uenced by it. Moreover, linear model without J 2 corrections is formulated and studied to see how the change in the reference orbital period and nodal drift errors aect the optimal solution. Finally, optimization results obtained with Hill's equations only are presented and compared. Figure 4.2: Optimal trajectory for the Full Linear Model (i ref = 70 ) 72 For the minimum-propellant-control rendezvous scenario set in Section [4] the pseudospectral method yields the solution shown in Figure [4.2]. Additional planar views of the optimal trajectory is shown is Figure [4.3]. Figure 4.3: Planar view of the optimal trajectory for Full Linear Model (i ref = 70 ) The green marker shows the initial position of the chaser with respect to the target. The nal maneuver time as discussed in the previous section is xed at one orbital period of the reference orbit which when calculated comes to t f = 1:72 hr and ends at the red marker. It is to be noted that the blue markers on the trajectory are the points where the thrust is applied and it is assumed that the thrust applied is impulsively. 73 Figure 4.4: Position and Velocity proles for the Full Linear Model (i ref = 70 ) Figure 4.5: Minimum propellant control proles for the Full Linear Model (i ref = 70 ) 74 Figure [4.5] shows the time history of the relative position and velocity of the optimal trajec- tory. Along with that, the control prole for the three normalized components of the thrust force is displayed in Figure [4.4]. It is clearly seen that as there is no signicant cross-track motion, the f z control is not very active. On the other hand because of the large amplitude of the motion in the radial direction (^ y) and the two rapid consequential change in the radial directions makes the f y control very active as it maxes out on the thrust bounds. It should be noted that the relative tolerance of 10 3 was set to solve the dierential equations. Addition to that the terminal tolerance of the order of 10 7 is set for the NLP solved by GPOPS, which does not necessarily indicate the quality of the solution. The reason is that for the solution the parameters of the trajectory are only being computed at 57 nodes only where the optimality conditions are enforced. For the implementation of the Gauss Pseudospectral method using GPOPS-II the guess for the initial states corresponds to the initial conditions and the guess for the nal state consisted of zeros for the entire state. The guess for the control history was zero at the initial and nal times for all controls. The solution obtained appears to be feasible but can be used only for oine computations, i.e. on an open-loop guidance scheme. 4.1.1 Eect of Non Corrected J 2 Perturbation Figure [4.6] adds to the study in this section by showing the eect of the dierential J 2 without corrections for 70 inclination on the optimal trajectory prole. This simulated orbit as a part of the optimal solution is alongside plotted with optimal trajectories of the Full Linear Model that implements all J 2 corrections and the Hill's equations so as to fully understand the sizable errors in the optimal positions caused by the the orbital period mismatch between the target and the chaser orbits and the nodal drift caused by the separation of the ascending node. 75 Figure 4.6: Optimal guidance trajectory for dierent models (i ref = 70 ) 4.1.2 Eect of Inclination on the Perturbed Optimal Trajectories As the eect of J 2 potential is directly proportional to the inclination of the target of the reference orbit. A generalized plot below in shown in Figure [4.7] that shows the J 2 perturbation eect on the optimal trajectories as the inclination of the reference orbit is increased. Figure 4.7: Optimal trajectory for the full linear model for dierent inclination It is seen that at lower inclinations the optimal route of the chaser is very similar to that of the Hill's equation which should be the case as we already discussed that the eect of J 2 is not 76 very signicant at inclinations near to equator (0 ). For all the simulation plotted in Figure [4.7], the same rendezvous scenario is considered with xed nal maneuver time of one orbital period. 4.2 Minimum Control Solution for Hill's EOMs This sections draws dierences in the optimality between the Hill's Equations and the Full Linear Model that integrates both J 2 with corrections and quadratic drag perturbations. The following Figure [4.10] shows the optimal trajectory for the Hill's equations. As expected a solution similar to Hohmann transfer is achieved. It can be hereby seen that the use of the Hill's equations for the study of optimal trajectories for rendezvous and proximity operations presents very conservative results. Figure 4.8: Optimal trajectory for the Hill's Equation and the Full Linear Model (i ref = 0 and = = 0). The above Figure [4.10] also shows the propellant-optimal trajectories of the full linear model for 0 inclination and assuming the drag constant of both the spacecrafts as zero i.e. = = 0 along with the the linearized Hill's equations of motion. Comparing the plotted optimal trajecto- ries validates the correctness of the formulated linear model relative equations. An insignicant 77 error noticed in the relative position is in the same order of magnitude as seen in Figure [3.23] corresponding to the non circular reference orbit used for the derivation of the full linear model. Figure 4.9: Minimum propellant control proles for the Full Linear Model for i ref = 0 and = = 0 Figure 4.10: Minimum propellant control proles for the Hill's Equation Figure [4.9] and Figure [4.10] shows a plot of the resulting control history of the three applied thrust forces corresponding to the optimal trajectory shown in Figure [4.8] using the linearized 78 Hill's equations of motion and the full linear model with perturbations at 0 inclination and zero drag constants. Comparing the identical thrust proles furthermore validates the correctness of the derived full linear model with perturbations for this study. 4.3 Cost Function Having established what the proper cost function is in the section above, a reasonable question to ask is the penalty incurred for using the wrong cost. Before discussing the penalty in fuel con- sumption it is quite important to note an important penalty in the design of the control system itself. The thrust force (or acceleration) appears linearly in a Newtonian dynamical system: this is a direct consequence of Newton's Laws of motion and not a simplication from linearization. In minimizing such control ane systems, barring the possibility of singular arcs, the L 1 -optimal controller has a bang-o-bang structure. That such controllers are indeed minimum-fuel con- trollers is well-known to astrodynamicists. On the other hand, quadratic-cost-optimal controllers are continuous con- trollers. Continuous controllers are frequently not desirable for spacecraft guidance and control since these controllers typically create undesirable eects on the payload. For ex- ample, thrusting increases the microgravity environment on the space station or induces undesirable eects on precision pointing payloads. Hence it is preferable to do much of the sci- ence during the \o periods". Thus, not only does the quadratic cost take more fuel, as will be shown presently, it creates new systems-engineering problems that were nonexistent prior to active control considerations. Rewriting the cost again for completeness, for a scaler valued function, f :R !R, the L p -norm of f, denoted byjjfjj p (<1) is dened by, jjfjj L p := Z jf(t)j p dt ! 1=p (4.7) 79 wherejj denotes the absolute value. For vector valued functions,f : R ! R n , n > 1, the L p -norm,jjfjj L p is frequently dened to be derived from Equation 4.7 withjj replaced by the Euclidean norm inR n . Thus, for example, if n = 2 so that f(t) = (f 1 (t);f 2 (t)), then, by this denition of a norm,jjfjj p L is given by, jjfjj L p := Z q f 2 1 (t) +f 2 2 (t) p dt ! 1=p (4.8) Applying this denition for our cost described in Section 4, Equation 4.4 J() J Q =jju()jj L 1 = Z t f t0 q u 2 x (t) +u 2 y (t) +u 2 z (t)dt (4.9) J Q =jju()jj L 2 = Z t f t0 (u 2 x (t) +u 2 y (t) +u 2 z (t))dt (4.10) Clearly Equation 4.10 described above, J Q =jju()jj 2 L 2 , the L 2 -norm of the cost does not measure fuel. On the other hand J Q =jju()jj L 1, the L 1 -norm of the cost does indeed measure fuel consumption. Figure 4.11: Optimal trajectory for the Full Linear Model (i ref = 70 ) 80 Figure 4.12: Minimum propellant control proles for the Full Linear Model for i ref = 70 with 6 eective controllers Figure 4.13: Minimum propellant control proles for the Full Linear Model for i ref = 70 with 3 eective controllers 81 As seen in Figure 4.11, the eect of the L 2 cost on the dynamics exaggerates the control eect required for the optimal trajectory. And the eect of number of controllers highlight the correctness of how transnational controllers should be implemented in the optimization problem. As any movement in one direction is facilitated by pair of controllers instead of a singular thrust control. It can be calculated from the impulse of the propellant used I sp = 70s along with the change in the V at each of the collocation points, fuel consumed for 6-controllers design for L 1 -norm is 1425.59g compared to L 2 -norm which is 2995.27g, which is an increase of +110:11%. 3-controller paints a similar story, where L 1 -norm is 1695.10g compared to L 2 -norm which is 1793.27g, which is an increase of +5:81%. Figure 4.14 highlights similar eect seen on the optimal solution for the Hills Relative motion model. The eect of controller choice is not signicant given the dynamics are accounting for any nonlinearities pertinent to J 2 and quadratic drag perturbations. Designing a controller under a circular earth assumption for 0 inclination orbits makes the choosing of cost norm and number of controller signicantly easier. Figure 4.14: Optimal trajectory for the Hill's Relative Motion Model (i ref = 0 ) 82 Figure 4.15: Minimum propellant control proles for the Hill's Model fori ref = 0 with 6 eective controllers Figure 4.16: Minimum propellant control proles for the Hill's Model fori ref = 0 with 3 eective controllers 83 There is a denite increase in the reported fuel usage under L 2 -norm, but not as material when accounting for motion under dierent perturbations. Fuel consumed for 6-controllers design for L 1 -norm is 329.23g compared to L 2 -norm which is 437.69g, which is an increase of +32:94%. Similarly for 3-controller design, L 1 -norm comes out to 371.26.g compared to L 2 -norm which is 390.47g, which is an increase of +5:17%. 4.4 Eect of Altitude and Maximum Allowable Thrust This section shows the eect of varying altitude on the optimal solution for dierent times-to- destination T f , that is the orbital period of the target. Figure 4.17 shows for smaller period sub 1T period the eect on the optimal solution and the fuel consumed is not that substantial. Although eect of drag quickly becomes signicant when the orbital periods are increased, which has been studied in Section 3. An interesting thing to note here is with increase in time-to-destination, as the optimization uses this xed time, chaser spacecraft could get into a perching state and keep expelling fuel until the desired time-to-destination constraint is met. If one would optimize also on the optimal tile while minimizing control cost, given the results, the solution would not be far from 1T period . Similar to altitude, as seen in Figure 4.18 the eect of maximum allowable thrust on the optimal trajectories for smaller orbits is not as signicant. This would make sense as the optimal solution for the lowest allowable thrust does not saturate itself, this can be seen from the control proles shown in Figure 4.4. Only one of 3 controllers saturate at 8:33 mm/s 2 forT max = 50mN and that too at the very beginning and during the end of the control prole, signifying that its primarily maxing out to give itself a head-start at the beginning and once more for stopping itself to reach the desired state of zero nal velocity. So any additional allowable trust would not be used and the nal solution would be fairly similar for smaller time-to-destination constraint. Although given this logic, for long periods one would then expect optimal solution using more 84 thrust at each point. This is clearly seen in the results, for orbital periods greater than 1, for higher maximum allowable trust, the chaser spacecraft is much close to the target than its lower trust counterpart. Figure 4.17: Minimum propellant trajectories for variable altitude and dierent time-to- destination constraint 85 Figure 4.18: Minimum propellant trajectories for variable altitude and dierent time-to- destination constraint 86 Chapter 5 Linear Quadratic Regulator This Chapter builds on previously studied minimum eort controller. LVLH (^ x^ y^ z) coordinate system described in Section [2.1.3] is used for the model formation. Similar to last chapter, we start from the arbitrary relative position and would like to bring the two spacecraft together for proximity operation. Given a controllable system, and either measurement of the full-state or an observable system with a full-state estimate, there are many choices of stabilizing control laws of form ~ u t . It is possible to make the eigenvalues of the close-loop system A BK arbitrary stable, placing them as far as desired in the left-half of the complex plane. However, overly stable eigenvalues may actually require exceedingly expensive control expenditure and might also result in the actuation signals that exceeded maximum allowable values. Choosing very stable eigenvalues may also cause the control system to over-react to noise and disturbances, much as a new driver will over-react to vibrations in steering wheel, causing the close loop system to jitter [nathon texbook]. Over stabilization can counterintuitively degrade robustness and may lead to instability if there are small time delays or remodeled dynamics. Choosing the best gain matrix K to stabilize the system without expending too much control eort is an import goal in optimal control. A balance must be struck between the stability of the 87 control-loop system and the aggressiveness of control. It is important to take control expenditure into account • To prevent the controller from over-reacting to high-frequency noise and disturbances • Actuation does not exceed maximum allowed amplitudes • Control is not prohibitively expensive. The above is encapsulated in dening a proper cost function that has the ability to penalize the state & control as desired. J(t) = Z t 0 x(t)Qx T (t) + u(t)Ru T (t)dt (5.1) The above cost function equation [5.1], balances the cost of eective regulation of the state with the cost of the control. The matrices Q and R weight the cost of deviations of the state from zero and the cost of actuation, respectively. The matrix Q is positive semi-denite, and R is positive denite; these matrices are often diagonal, and the diagonal elements may be tuned to change the relative importance of the control objective. Adding such a cost function makes choosing the control law a well-posed optimization problem. The linear quadratic regulator (LQR) control law u =K r x is designed to minimize J(t) as time t f goes to ininity, J = lim t!1 J(t). LQR is so-named because it is linear control law, designed for a linear system, minimizing a quadratic cost function, that regulate the state of the system to lim t!1 x(t) = 0. Because the cost-function is quadratic, there is an analytical solution for the optimal controller gains K r , given by K r = R 1 B X (5.2) where X is the solution to an algebraic Riccati equation: 88 A X + XA XBR 1 B X + Q = 0 (5.3) The LQR controller is shown schematically in Figure [5.1] below Figure 5.1: Schematic of the liner quadratic regulator (LQR) for optimal full-state feedback. The optimal controller for a linear system given measurements of the full state, y = x is given by proportional control u =K r x where K r is optimal gain at a given time t 5.0.1 Choosing Q and R matrices Choosing proper Q and R is important for optimization convergence. The following implemen- tation was used as a baseline for minimizing the cost J(t). As seen in equation [5.4] and [5.5] the diagonal values for both Q and R are scaled with a factor of their inverse maximum allowable constraint. This implementation is a classic Bellman's implementation, it accounts for applying "same" magnitude of penalty to each state relative to their maximum allowed values. For dierent units along the diagonal this makes the cost penalty unit less. 89 Q = 2 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 4 Q 1 x 2 max : : : : : : Q 2 y 2 max : : : : : : Q 3 x 3 max : : : : : : Q 4 _ x 2 max : : : : : : Q 5 _ y 2 max : : : : : : Q 6 _ z 2 max 3 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 5 (5.4) where, Qi are the scaling parameters. This can be used to apply a more robust cost penalty on a specic state variable. Like Q matrix, R also uses a similar approach, Ri here acts as a scaling factor that let us penalize cost on each controller. These matrices can be changed dynamically by using a time varying value for both Qi and Ri . R = 2 6 6 6 6 6 6 4 R 1 u 2 xmax : : : R 2 u 2 ymax : : : R 3 u 3 zmax 3 7 7 7 7 7 7 5 (5.5) 5.0.2 Numerical Solutions to the LQR problem There are many possible solutions for solving of LQR cost as the problem truly depends on the chosen Q and R. Choosing a Q matrix that is time varying versus one that is not is a controls engineers design choice. As discussed in the previous section penalty matrices imple- mented are normalizing each state at any given time by its maximum allowable vales; state bounds/constraints. For simplicity, given the problem is already complex, no scaling parameters are chosen to individually actuate cost on the each state. The nal form of Q matrix is shown in Equation [5.6] 90 Q = 2 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 4 1 x 2 max 0 0 0 0 0 0 1 y 2 max 0 0 0 0 0 0 1 x 3 max 0 0 0 0 0 0 1 _ x 2 max 0 0 0 0 0 0 1 _ y 2 max 0 0 0 0 0 0 1 _ z 2 max 3 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 5 (5.6) Resonable values for state bounds where chosen from the optimal fuel controller solution that were implemented in the previous section. From the position and velocity proles as seen in Figure [4.5], j~ xj =)jxj;jyj;jzj 5 km Position Bounds j _ ~ xj =)j _ xj;j _ yj;j _ zj 8 m/s Velocity Bounds For the R matrix, scaling factor was chosen as Ri =r g . Here r g is the range-to-goal factor. For the optimization problem described here this is simply the relative distance of the chaser from target at any given point in time. The denation of the R matrix along withr g is shown in Equation [5.7]. Also to note, the penelty also accounts for the maximum allowable trust or control eort in each direction for their respective thrusters. This makes the penelty on the control time varying and dynamic, i.e. as the chaser gets closer to the targer, the eect on the use of the thruster decreases. It is when the chaser is furthest from the targer, most control/thrust is expelled. 91 R = 2 6 6 6 6 6 6 4 rg u 2 xmax 0 0 0 rg u 2 ymax 0 0 0 rg u 3 zmax 3 7 7 7 7 7 7 5 (5.7) where r g = p x 2 +y 2 +z 2 Maximum bounds for the control eeort by the thrusters is already explained in previous section. For the mission design discribed by U.Bogdan et. al. in [31] it is assumed that each thruster produces 50 mN, with the chaser mass as 24 kg - control bounds can be calculated by F T =m chaser ju i j =)ju x j;ju y j;ju z j 8:33 mm/s 2 Control Bounds 5.1 LQR Control Solution This section highlights control solutions for the implemented Q and R matries above. Figure [5.2] shows the LQR optimal trajectory for 1-orbital period. As seen in the gure, LQR solution closely resembles for optimal fuel solution. Given LQR is a close-loop controller, this fact is highlighed in the position and velocity proles of the optimal solution as seen in Figure [5.3]. Close loop controller tries to minimize the error but damping the states to desired values as quickly as possible similar to a close-loop step respose solution. 92 Figure 5.2: LQR Optimal trajectory for CHW Model (1 orbital period) Figure 5.3: Position and Velocity proles for CHW Model (1 orbital period) 93 Thrust prole as seen in Figure [5.4] validate the chosen R matrix. As ther g value gets smaller, i.e. the relative distance between the chaser and the targer gets smaller, less fuel is explled. Figure 5.4: LQR optimal trajectory control prole for CHW Model (1 orbital period) 5.1.1 Variable t f This section implements varibale time for the optimization routine. It was seen from numerous simulations that minimizing LQR cost for a xed window of simluation period (1 orbital period) the optimization problem becomes more strict and constrainted. The solution itself with xed simulation period becomes very sensitive to the chosen maximum bounds for the state variables and control eort. It became apparent that easing this time constraint would make the solution more robust, and get the optimization convergence to better tolerances and faster. The solution below uses the same Q and R matrices as implemented above. The only dierence being unbounded simulation time. The control solution now will converge at t f when the cost is minimized with higher tolerances. 94 Figure 5.5: LQR Optimal trajectory for CHW Model (2.572 orbital period) Figure 5.6: Position and Velocity proles for CHW Model (2.572 orbital period) As seen in Figure [5.5] the solution looks much smoother compared to the solution that was constrained to 1 orbital period as seen in Figure [5.2]. The optimization problem converges at t f = 2:572T 0 95 The position and velocity proles showen in Figure [5.6] validate the close-loop nature of LQR cost. 5.1.2 Comparing v This sections compares the v for both minimum control and LQR solution. Following summary shows v for dierent senarios modelled. It is clearly seen that v LQR comes ahead of the optimal fuel controller v MC . This is the penelty of implementing a more design friendly controller. v MC v LQR 1T 0 9.36 m/s 13.23 m/s 2:572T 0 12.21 m/s 14.10 m/s In order to control penelty on state variable, which would be necessary a requirement for path planning. The added tuning of the Q and R matrix makes the LQR solution not the most fuel ecient. 96 Chapter 6 Neuromorphic Controller Computation cannot be physically realized without time. Whether we consider the discrete switch- ing of voltages across transistors in a digital computer, or the continuous ow of charges across ion channels within the brain { all physical systems, that we currently know of, make progress on their computations by changing over time. Digital computers pipeline their computations through multiple cores with energetically-expensive access to distant memory, while our brains consume only twenty watts of power by harnessing the dynamics of a biophysical substrate.A leitmotif of this section is a computational paradigm in between these two extremes: neuromorphic comput- ing. The goal is to draw inspiration from how the brain dynamically performs computations to engineer low-power digital and analog circuits that emulate its fundamental principles. Despite growing excitement, there exist many challenges when it comes to using neuromorphic hardware to perform some desired task. Generally speaking, neuromorphic programming forces us to rethink our traditional, Von Neumann (1958), understanding of computation, and redesign our algorithms to leverage distributed networks of noisy, heterogeneous units, that locally process their inputs and communicate via spikes { also known as Spiking Neural Networks (SNNs). These networks serve as powerful models of neural computation, capable in theory of everything that we can do, while consuming extraordinarily low amounts of energy. Biological nervous systems already leverage the computational power of SNNs quite successfully. Thus, our approach is to 97 borrow existing principles and models from neuroscience, exploit their properties using tools from mathematics and engineering, and use software methods to synthesize the resulting networks in silico. Moreover, the frameworks and methods that we use to accomplish this should be extensible to incorporate increasingly-sophisticated levels of biological detail so long as such details are computationally useful at a functional level, or shown to benet some resource-power-accuracy trade-o. By virtue of emulating core principles of brain function, we hope to realize algorithms that consume far less power, and scale more readily to massively parallelized computing architectures running in real-time. But such algorithms must also be practically useful. To ip this on its head: we require frameworks that describe what some hardware is doing at the lowest level, so we may systematically harness its capabilities to perform useful computations at the highest level. Anal- ogous to how a single compiler can translate the same C++ code to dozens of physically distinct digital architectures, we strive to systematically describe how the same dynamical computations can be exibly mapped onto distinct spiking neuromorphic architectures. The general approach addressed in this section adopts a software stack and nonlinear dynamical systems framework for describing computations within vector spaces. This approach is known as the Neural Engineering Framework (NEF; Eliasmith and Anderson, 2003). The general-purpose neural simulator, Nengo (Bekolay et al., 2014), natively supports NEF networks and a variety of neuromorphic backends. The ultimate goal of this approach is to develop a working model of simpler dynamics that support and pave the way into implementing neuromorphoic controllers for space systems. Although we still have a long ways to go, this section aims to highlight a promising path for neuromorphics. When confronted with the ambitious goal of articially reproducing human intelligence, com- puter scientists will often point towards the popularity of deep learning (LeCun et al., 2015) and its recent successes such as AlphaGo (Gibney, 2016). First and foremost, we seek to embrace the methods of deep learning, by incorporating its architectures and training methodologies into our 98 toolbox. This is in fact the mandate of Nengo DL (Rasmussen, 2018), a Python package that combines Nengo and TensorFlow into a hybridized framework that auto- mates backpropagation across spiking, non-spiking, and NEF-style recurrent neural networks broadly referred to as the class of Articial Neural Networks (ANNs). Deep learning methods have found enormous success, due in part to the black-box nature of the training procedure, which requires a programmer to have virtually no knowledge of the underlying problem being solved. However, we do not believe that deeper networks, better training heuristics, nor larger datasets, will be enough to achieve our aforementioned goal. Some concerns for conven- tional deep learning approaches include: the amount of power required to scale traditional ANNs to the size of the human brain (Furber, 2012), the diculties tuning hyperparameters at scale (Bergstra et al., 2015), and the abundance of adversarial attacks on networks trained via backpropagation (Su et al., 2019). More generally, the \no free lunch theorem" informs us that a black-box learning algorithm, on its own, will never be enough (Wolpert, 1996); inevitably, it becomes necessary to incorporate a priori knowledge into the structure of ANNs in a manner that biases the solution space of the training procedure towards the problems being solved. But even more fundamentally, the majority of ANNs in use today, including AlphaGo, are applied to static inputs to produce static outputs. For a board game such as Go, this is forgiv- able, as the entire state of the game is encapsulated by the current state of the board itself, and so there is technically no need to consider any history in order to optimize across future states. Nevertheless, this points towards a thematic concern: time is not being treated as a continuous dimension of computation in many of the most successful methods, which may limit their applica- bility to many real-world problems. For instance, convolutional ANNs are most often deployed in video processing tasks by classifying objects in single frames, sampled at some discrete interval, and independently from one another. Recurrently connected networks are an exception that we embrace, but they come with many challenges. In general, most methods do not respect time as 99 being a physical dimension that underpins the dynamics of the system's input, output, or pro- cessing units; in typical ANNs, time is not a variable in the neuron models, synapse models, nor any of the computational elements. In contrast, the human brain is a physical system that continuously processes information over time. This is crucial for a cognitive agent embodied in a dynamic environment wherein all interactions depend on internally and externally evolving states. The brain naturally learns from interacting with its environment over many dierent time-scales. For instance, our brains are not programmed to recognize static objects at discrete moments in time. Rather, we have evolved to interact with objects in continuously changing environments, from which the perception of object recognition emerges. From a very young age, children can intuitively infer the physics of moving objects, and learn how to interact with those objects to achieve desired outcomes. This continues through adulthood, up to the level of socially coordinating our thoughts and behaviours, sculpted by our cumulative experience with the world. This occurs so naturally that we often take it for granted, but at a very basic level, our brains are necessarily constructing explicit and implicit dynamical models of reality, and applying these models to perform behaviourally relevant computations in real-time. Consider a real-world example such as driving a car. We must constantly integrate new information with our existing understanding of where we are, where we are headed, and what might happen along the way. This requires not only an intuitive understanding of the physics of the car, but models of how other drivers behave on the road in the context of changing trac lights and road conditions, all of which play out in parallel while we exibly plan and coordinate our actions. Likewise, consider any actively engaging task or hobby that you enjoy. I speculate this activity requires some degree of mentally or physically coordinated interaction with the world, on a similar level of dynamic complexity as in the driving example, such as: playing a video game, participating in a sport, engaging with some form of media, playing music, drawing, dancing, and 100 so forth. Besides such activities being meaningful, fun, and beautiful, they all share the common motif of recruiting a variety of systems in a dynamically coordinated fashion. To help appreciate how extraordinary such tasks are for the human brain, we should recognize that the mechanisms involved do not have a perfect memory of past inputs. Each neuron locally processes some signal by continuously encoding it into spike trains. These action potentials evoke postsynaptic currents (PSCs) that decay exponentially on the order of milliseconds. Memories emerge from these dynamics, as changes to activation patterns, adaptive states, and synaptic weights, all re ect the histories of representations along a myriad of time-scales. The time-scales over which these mechanisms interface with the world, and interact with each other, determines, and constrains, all of our perceptions, thoughts, and actions. Remarkably, this is all accomplished by a collection of approximately one-hundred billion neurons, connected by over one-hundred trillion synapses, processing these signals sparsely across space and time, while consuming only twenty watts (Koch, 2004) { about the same energy as a typical household CFL bulb. This "signal processing view" of computation is fundamentally dierent from conventional views of computation targeted for conventional computers. In the former view, time is intrinsic to the state of all physical "devices" performing the computation (i.e., neurons and synapses), which are themselves necessarily tied to the timing and dynamics of the overall function as a consequence of actually implementing the function themselves. However, in the conventional view, the time that a computation takes is an implementation detail that is mostly decoupled from the input- output relation. As such, there is currently an important distinction between the computational approaches taken by the neuromorphic community and by those who program traditional digital computers. This is not simply a matter of biological detail, but more crucially a matter of thinking about the system that is actually solving these tasks as a collection of signal processors whose dynamics are coupled to the task at hand by the physics of time. If our aim is to build articial systems that are as useful, robust, capable, and scalable as the human brain, then we should strive to understand computation in such terms. 101 6.0.1 Neural Engineering Framework NEF provides a quantiable way to represent neurons. It explains how neurobiological systems represent the world, and how they use those representations, via transformations, to guide be- havior. In order to precisely dene the representation relation we turn to the quantitative tools of engineering. Specically, we believe that there is a close tie between neural representations as understood by neuroscientists and codes as understood by communications engineers (Reza 1961). Codes, in engineering, are dened in terms of a complimentary encoding and decoding procedure between two alphabets. For example, Morse code is dened by the one-to-one relation between letters of the Roman alphabet, and the alphabet composed of a standard set of dashes and dots. The encoding procedure is the mapping from the Roman alphabet to the Morse code alphabet and the decoding procedure is its inverse (i.e., the mapping from the Morse code alphabet to the Roman alphabet). Biological brains possess a signicant amount of structure, with neuroanatomical constraints that dier considerably between functional areas. We interpret this observation pragmatically, as a call to: (1) incorporate prior knowledge of a desired set of tasks into the network's construc- tion, and (2) constrain a network's function using available models of the underlying biophysical mechanisms. However, the standard formulation of the NEF makes a number of simplifying as- sumptions regarding these mechanisms. For example, the NEF typically assumes an exponential model of the postsynaptic current (PSC) evoked by an action potential, which has a biologically- implausible, instantaneous, rise-time. This model is also known as a rst- order lowpass lter. In contrast, the synapse models in mixed-analog-digital neuromorphic chips are non-ideal, featuring higher-order dynamics due to parasitic capacitances (Voelker et al., 2017a). Similarly, the synapse models commonly used in biological models incorporate distinct rise-and fall-times due to sepa- rate time-scales of transmitter binding and unbinding, as well as axonal transmission delays due to the nite-velocity propagation of action potentials (Roth and van Rossum, 2009). In digital 102 hardware, spike-transmission delays are congurable (Davies et al., 2018) and necessary at the time-scale of the simulation timestep (Bekolay et al., 2014). In the context of RNNs, the NEF provides generic machinery for training the weights of a recurrent network of spiking neurons in order to implement some particular dynamical system (Voelker and Eliasmith, 2016). This is made ecient through three important steps: (1) weight matrices are factored into nq encoding and decoding matrices, where q2 (1) is the dimen- sionality of the dynamical state-vector, and n is the number of neurons, (2) the optimization problem is recast as a batched least-squares problem that may be solved oine without needing to explicitly simulate the network over time, and (3) the dynamics of the synapses are leveraged as temporal bases for the overall system. Conceptually, the NEF consists of three mathematical principles used to describe neural com- putation: (1) representation, (2) transformation, and (3) dynamics. The NEF is most commonly applied to building dynamic (i.e., recurrent) spiking neural networks, but also applies to non- spiking and feed-forward networks. Its primary strengths reside in providing a well-understood and ecient means of engineering spiking neural models, and programming neuromorphic hard- ware, to approximate mathematically described computations with precision that scales as ( p n) in the spiking case and (n) in the non-spiking case (Eliasmith and Anderson, 2003; Boahen, 2017). We now provide an overview of these methods, applied to training both feed-forward and recurrent connection weights, with an emphasis on mapping linear dynamical systems onto SNNs. 6.1 Principle 1 { Representation Let x(t)2R q denote aq-dimensional continuous-time varying signal, that is to be represented by a population ofn spiking neurons. To describe this representation, we dene a nonlinear encoding and a linear decoding that together determine how neural activity relates to the represented vector. 103 NEF implements encoders E = [e 1 ;:::; e n ] T 2R nq , gains i > 0 and biases i(i = 1:::n) as parameters for the following encoding s i (t) = i (e i x(t)) + i a ( t) =G i [s i (t)] = m (tt i;m ) (6.1) wheres i (t) is the input current to thei th neuron,a i (t) is the neural activity generated by the i th neuron encoding the vectorx(t) at time t,G i [] is the nonlinear transfer function modelling the spike-generation of a single neuron, () is the Dirac delta, and t i;m is the sequence of spike-times generated by the neuron model in response to the input current. A common choice of neuron model for G i [] is the leaky integrate-and-re (LIF) neuron model (Koch and Segev, 1998). This model maintains a one-dimensional voltage trace,v i (t), by continu- ously integrating a current-sources i (t) over time, with some leakage governed by the time-constant RC { functionally equivalent to a resistor-capacitor (RC) circuit, a rst-order lowpass lter, or a convolution with an exponentially decaying impulse response. An action-potential, or \spike," is generated when the voltage reaches a threshold, in which event the voltage is reset to its resting potential, where it remains clamped for the duration of some refractory period, ref . The heterogeneous parameters that dene the encoding, (e i (t); u ; i ) : i = 1::n determine the variety of nonlinear responses from the population. Typically the length of each encoder, jje i jj 2 , is xed to 1, while i > 0 controls its length. Then e i may be interpreted geometrically as a \preferred direction" vector, such that its dot product similarity with x(t) determines the relative magnitude of s i (t).The bias terms i eectively determine the sparsity of representation (i.e., which neurons are active for a given x(t)), while the gain terms i determine the density of spikes relative to ref (i.e., how many spikes are generated by an active neuron). These parameters 104 can be randomly sampled from distributions constrained by the domain of x(t) and the dynamic range ofG i [], t from neuroanatomical data (e.g., known tuning curves, preferred directions, ring rates, sparsity, etc.; see for instance Friedl et al. (2016)), predetermined using prior knowledge of the desired transformation (Gosmann, 2015), or trained via backpropagation through time (Rasmussen, 2018). By default, we select these parameters such that the neural responses are uniformly distributed with respect to x across the q-dimensional hyperball (Gosmann, 2018, pp. 51{65). In essence, equation 6.1 denes a high-dimensional nonlinear projection of the vector x(t), by taking its dot product with an encoding matrixE, and injecting the result inton heterogeneous spike-generators. Figure [6.1] shows the lifecycle of the input signal, undergoing the encoder (turing curves) which at threshold create voltage spikes. Figure 6.1: Principle 1 & 2 of NEF in action. The life-cycle of the input signal how its encoded and decoded back to originality 105 6.2 Principle 2 { Transformation The second principle of the NEF addresses the issue of computing transformations of the repre- sented signal. The encoding remains dened by equation 6.1 However, we now decode a ltered version of some function, f : R q R q , by applying an alternate matrix of decoders to the same ltered spike trains: D f = h d f 1 ;::::; d f n i f(x?h)(t)(a i ?h)(t)d f i = n i=1 m h(tt i;m )d f i (6.2) Now we must solve for D f . For both Principles 1 and 2, we optimize for D f over the domain of the signal, S = x(t) :t 0, which is typically the unit q-ball v2R q :jjvjj 2 1 or the unit q-cube [1; 1] q . This can be done, as is, by simulating the population over time and solving a least-squares problem. D f =argmin D2R nq Z S f(v) m i=1 (r i (v) +)d i 2 2 d q v (6.3) Note that this optimization only depends on r i (v) for v2 S, as opposed to depending on the signal x(t). In other words, the optimization is determined strictly by the distribution of the signal, and not according to its particular dynamics. Since we have replaced the optimization problem involving spikes (equation 6.2) with one involving static nonlinearities, this sidesteps the need to explicitly simulate the neurons over time when training. Furthermore, this is a convex optimization problem, which may be solved by uniformly sampling S (Voelker et al., 2017b) and applying a standard regularized least-squares solver to the sampled data (Bekolay et al., 2014). In the context of machine learning, the samples of S are the example inputs from training data. Similarly in the context of space dynamics, S become the desired/nal state conditions. 106 The accuracy of this approach depends on the quality of the following approximation: x(t) = v =) m i=1 (a i ?h)(t)d f i n i=1 ((r i (v +)?h)(t)d f i (6.4) In plain words, this presupposes that substituting a ltered spike-train with its correspondingly ltered (but noisy) rates does not change the optimal set of decoders for a target set of postsynaptic currents. A subtle, yet critical, observation is that we do not require the spike rates of individual neurons to re ect their idealized instantaneous rates. Rather, we exploit the fact that ltering and weighting spikes across the entire population will satisfy such an approximation (Eliasmith and Anderson, 2003, pp. 132{136). Implementing such a encoding and decoding facilitates approximations of any function. As n number of neurons increases i.e. n!1 the approximation of the input signal v(t) theoretically should be exact. Figure 6.2: Approximation of dierent input signals 107 Figure 6.2 shows reimplements the ramp-up input signal as shown in Figure 6.1 but with more neurons in the ensemble. It is seen clearly that the eect of the noise on the approximated decoded signal is signicantly lower when number of neurons are increased. This idea of transformation can be extended by applying another function to the decoded output. Figure 6.3 below shows a crude working of a decoded that is approximating the "square" of the input signal, albeit the approximation is not very exact but this can be resolved by simply increasing the number of neurons in the ensemble population. Figure 6.3: Approximation of dierent input signals 6.3 Principle 3 { Dynamics Principle 3 is a method of harnessing the dynamics of the synapse model for network-level infor- mation processing. We begin by focusing our discussion of NEF dynamics on the neural imple- mentation of continuous, linear time-invariant (LTI) systems _ x(t) =Ax(t) +Bu(t) y(t) =Cx(t) +Du(t) (6.5) 108 where the time-varying signal x(t) represents the system state, _ x(t) its time-derivative, y(t) the output, u(t) the input, and the time-invariant matrices (A;B;C;D) fully describe the system (Brogan, 1991). This form of an LTI system is commonly referred to as the state space model, but there are many other equivalent forms. Figure 6.4: Block diagram for an LTI system. The integrator is driven by the signal _ x(t). Repro- duced from Voelker and Eliasmith Figure 6.5: Block diagram for an LTI system, with the integrator replaced by a rst-order lowpass lter. Reproduced from Voelker and Eliasmith For LTI systems, the dynamical primitive that is, the source of the dynamics is the integrator as shown in Figure 6.4. However, in the model that we consider under the NEF described in 109 Figure 6.5, the dynamical primitive at our disposal is the leaky integrator, given by the canonical rst-order lowpass lter modeling the synapse which looks like , h(t) = 1 e t =L 1 n 1 s + 1 o (6.6) whereL 1 denotes the inverse Laplace transform, and the latter representation is referred to as a transfer function. Our approach is to represent the stat vector x(t) in a population of spiking neurons (Principle 1), such that this same vector is obtained by ltering an appropriate transformation with a leaky integrator (Principle 2). Thus, the goal of Principle 3 is to determine the transformations required to implement Equation 6.5, given that the required x(t) is obtained by some convolution with a leaky integrator, rather than the perfect integrator depicted in Figure 6.4. Principle 3 states that in order to implement equation 6.5 in a population of neurons that represents x(t), we must compensate for the eect of \replacing" the integrator with a leaky integrator (compare Figures 6.4 and 6.5) by driving the synapse with _ x(t) + x(t) instead of only _ x(t) This compensation is achieved as follows: implement the recurrent transformation (A + I)x(t), and the input transformationBu(t), but use the same output transformationCx(t), and the same pass through transformation Du(t) (Eliasmith and Anderson, 2003). Specically, this may be implemented in a spiking dynamical network by representing x(t) via Principle 1 and then using Principle 2 to decode the needed transformations. The resulting model is summarized in Figure 6.5. Notably, both systems compute the exact same signals x(t) and y(t); the dierence is the dynamical primitive, and the transformations needed to account for this change. The neural architecture that realizes this model consists of a population of neurons representing x(t), recurrently coupled to itself, with weights implementing the above transformations, and each synapse ltering its spike trains via convolution. The outcome 110 is that the desired state vector, x(t), is linearly projected onto the postsynaptic current of each neuron. Principle 3 is useful for accurately implementing a wide class of dynamical systems - integra- tor, oscillators, attractor networks, controllers, and even continuous memories to solve specic problems that frequently arise in neural modeling (e.g., Eliasmith and Anderson, 2000; Singh and Eliasmith, 2004; Eliasmith, 2005; Singh and Eliasmith, 2006). Using the principles outlined above, the following gure 6.6 shows the implemented inverted pendulum system that is using a neuromorphic adaptive PI controller. Figure 6.6: Adaptive PI controller for Inverted Pendulum system in Nengo GUI 111 Chapter 7 Conclusions 7.1 Summary The popularity of rendezvous and docking missions is growing on several fronts. Future missions including satellite servicing and repair, orbital debris removal, and on-orbit assembly all require rendezvous and docking. Several of these missions also require maneuvering with the target of interest after docking. Majority of the future missions are going to be relying on small satel- lites given the advent of many launch vehicles facilitating passenger payloads for reduced cost. Considering the size of these small sats, it becomes imperative to have an accurate description of the dynamics aecting the motion of the spacecraft. This thesis has focused on this crucial topic, a more higher delity linear in nature relative motion equations are developed. A generic control problem formulation is created to dene the purpose and application of the developed linear model. Today's rendezvous and docking missions are planned to very well-known, designed- to-be- docked-to targets. With these new mission types on the horizon, there also comes a higher degree of uncertainty in the system. When performing orbital debris removal or servicing of broken satellites, the state of the target object may be unknown. Tracking from ground sensors may not be able to resolve this uncertainty, and thus, missions must be able to perform well while 112 considering the lack of concrete knowledge of system parameters. In the context of rendezvous, docking and joint maneuvering missions, this thesis studies the eect of these uncertainties, its magnitude for various design parameters. Comparisons are drawn to the current state-of-the art models that are used for controller designs. Eect of these nonlinear perturbing forces on the optimal solution is also analyzed highlighting the need for their implementation when designing missions or controllers for rendezvous and docking missions. 7.2 Contributions The optimal close range rendezvous problem for an ADR spacecraft to approach a tumbling upper stage is formulated and analyzed in this research. The following list highlights the contributions made part of this thesis. • The eect on the satellite motion due to J 2 perturbation is well modeled and analyzed. It was shown that a general time average of the J 2 does not account for all the nonlinearities, that additional corrections terms have to be incorporated to fully capture the nonlinear J 2 dynamics. • In the context of rendezvous and docking missions, especially when facilitated by small cube sats. It was shown that the in uence of quadratic drag having not modeled into dynamics can have serious implications. Perturbations dominate in the in-track motion of the satellite and can deviate the satellite orbit from its nominal counterpart only after a few orbital period. • Upon highlighting the above eects and the neediness for their implementation. A higher delity linearized equations of relative motion were developed. Time averagingJ 2 along with quadratic drag and other corrective terms are theorized to fully capture all the nonlinearities due to these perturbations. 113 • An optimal control problem that minimized propellant usage is theorized as an obvious application to the developed higher delity linear model. Various sensitivity studies were conducted, ranging for choosing the correct cost function, to eects of time-to-destination constraints along with other altitude, inclination and max thrust constraints. • Given the later focused more on open loop problem. A close-loop controller application was also designed on the linear models. LQR controllers were studied with varying cost penalty to state and control errors. v was compared to minimum-control solutions to understand the cost of manually tuning and controlling penalty for designing of the trajectories compared its open loop counterpart. • A formal primer on neuromorphic controllers was summarized. The underlying theory of Neural Engineering Framework was utilized to implement a working neuromorphic adaptive PI controller for stabilizing an inverted pendulum system. This research was a stepping stone for showcasing that such biological plausible controllers that are very power and memory ecient can be used as valid controllers for physical systems and even for the navigation and control of the next generation of autonomous spacecraft. 7.3 Future Work Further study is required which can incorporate the the full 6dof model of the chaser spacecraft, where both the rotational and transnational dynamics will be considered. Moreover, to make the capture scenario more realistic collision-avoidance condition can be imposed in form of a path constraint. Along with that, docking/capture-enabling conditions can also be incorporated by matching the chaser's and target's docking station position and velocity vectors along with other attitude constraints. Optimal time and optimal energy conditions can also be studies for the rendezvous scenario to truly nd the optimal trajectory for target capture. To study the sensitivity of the initial 114 conditions, control problem can be formulated analytically using an indirect method to better understand the numerical results. Other control application ranging from PID to model predictive controllers can be studied under correct dynamics. 115 Reference List [1] Aaron,R. "Dynamical Systems in Spiking Neuromorphic Hardware" University of Waterloo, Canada 2019 [2] Alfriend K.T., Schaub H., Gim D., "Gravitational Perturbations, Nonlinearity and Circu- lar Orbit Assumption Eects on Formation Flying Control Strategies", AAS Guidance and Control Conference, Breckenridge, CO, Feb. 2-6, 2000 [3] Bate R.R., Mueller D.D. and White J.E., Fundamentals of Astrodynamics, 2nd Edition, Dover Publications Inc., 1971 [4] Burak, O. "Model predictive control for spacecraft rendezvous and docking with uncooperative targets" Nanyang Technological University, Singapore 2020 [5] Curtis H., "Orbital Mechanics for Engineering Students", 3rd Edition, Elsevier, 2014 [6] David, L,. "Military Micro-Sat Explores Space Inspection, Servicing Technologies" Space.com, www.space.com/businesstechnology/050722 Xss-11 test html, Technology, 22 July 2005 [7] Dornheim, M. A., "Orbital Express to Test Full Autonomy for On-Orbit Service," Aviation Week & Space Technology, www.aviationnow.com/avnow/ news/channel awst story.jsp?id=news aw0605 06p1.xml, 4 June 2006 [retrieved 20 July 2009]. [8] Eliasmith, C. and Anderson, H. "Neural Engineering - Computation, Representation, and Dynamics in Neurobiological Systems" 116 [9] Friend, R. B., "Orbital Express Program Summary and Mission Overview," Proceedings of SPIE: The International Society for Optical Engineering, Vol. 6958, 2008, Paper 695803. doi:10.1117/12.783792 [10] Gaias G. and D'Amico S., "Impulsive Maneuvers for Formation Reconguration Using Rel- ative Orbital Elements", Journal of Guidance, Control and Dynamics, Vol. 38, No. 6, June 2015 [11] G. Boyarko, "Spacecraft guidance strategies for proximity maneuvering and close approach with a tumbling object," PhD, 2010. [12] G. Boyarko, Y. Oleg, and R. Marcello, "Real-time 6dof guidance for of spacecraft proximity maneuver- ing and close approach with a tumbling object," in 2010 AIAA/AAS Astrodynamics Specialist Conference, ser. Guidance, Navigation, and Control and Co-located Conferences. American Institute of Aeronautics and Astronautics, 2010. [13] G. Boyarko, O. Yakimenko, and M. Romano, "Optimal rendezvous trajectories of a controlled spacecraft and a tumbling object," Journal of Guidance, Control, and Dynamics, vol. 34, no. 4, pp. 1239{1252, 2011. [14] Humi M. and Carter T., "Fuel-Optimal Rendezvous in a Central Force Field with Linear Drag", AAS/AIAA Space Flight Mechanics Meeting, Santa Barbara, California, 2001 [15] Humi M. and Carter T., "Rendezvous Equations in a Central-Force Field with Linear Drag", Journal of Guidance, Contol and Dynamics, Vol. 25, No. 1, Jan.-Feb. 2002 [16] Humi M. and Carter T., "Clohessy { Wiltshire Equations Modi eed to Include Quadratic Drag",Journal of Guidance, Contol and Dynamics, Vol. 25, No. 6, Nov.-Dec. 2002 [17] Humi M., "Low-Eccentricity Elliptic Orbits in a Central Force Field with Drag", Journal of Guidance, Contol and Dynamics, Vol. 33, No. 5, Sept.-Oct. 2010 117 [18] Ichimura, Y., and Ichikawa, A., \Optimal Impulsive Relative Orbit Transfer Along a Circular Orbit," Journal of Guidance, Control, and Dynamics, Vol. 31, No. 4, 2008, pp. 1014{1027 [19] Lawdin D., \Optimal Trajectories for Space Navigation," Butterworths, London, 1963. [20] Michael, R. "How to Find Minimum-Fuel Controllers" AIAA Guidance, Navigation and Control Conference and Exhibit, Providence, Rhode Island 2004 [21] Michael, R. "A Roadmap for Optimal Control: The Right Way to Commute" Annals of the New York Academy of Science, January 2006 [22] Michael, R. "Space Trajectory Optimization and L 1 -Optimal Control Problems", Modern Astrodynamics, Elsevier Astrodynamics Series 2006 [23] NASA, DART Mission, http://www.nasa.gov/missions/ -science/dart into space.html/ [24] Prussing, J. E., and Conway, B. A., Orbital Mechanics, Oxford Univ. Press, New York, 1993, pp. 139{169. [25] Rao, A. V., Benson, D. A., Darby, C. L., Patterson, M.A., Franscolin, C., Sanders, I., and Huntington, G. T. (2008), "GPOPS: A MATLAB Software for Solving Multiple-Phase Optimal Control Problems Using The Gauss Pseudospectral Method" ACM Transsaction on Mathematical Software [26] Roscoe C.W.T., Westphal J.J., Griesbach J.D., and Schaub H., "Formation Establishment and Reconguration Using Dierential Elements in J 2 -Perturbed Orbits", IEEE Aerospace Conference, Big Sky, Montana, March 1-8, 2014 [27] Ross Michale I., "How to Find Minimum-Fuel Controllers", AIAA Guidance Navigation and Control Conference, 16-19 August 2004, Providence, Rhode Island [28] Ross Michale I., "Space Trajectory Optimization andL 1 -Optimal Optimal Control Prob- lems", Modern Astrodynamics, pp. 155-188, Elsevier Asytodynamics Series, 2006 118 [29] Sakawa Y. and Shindo Y., "Optimal Control of Container Cranes", Automotica, Vol. 18, No. 3, pp. 257-266, 1992 [30] Schaub H., Vadali S.R., Junkins J.L. and Alfriend K.T., "Spacecraft Formation Flying Con- trol using Mean Orbit Elements", Journal of Astronautical Sciences, Vol. 48, No. 1, Jan.-March 2000, Pages 69-87 [31] Schaub, H., and Alfriend, K. T., \J 2 Invariant Reference Orbits for Spacecraft Formations," Celestial Mechanics and Dynamical Astronomy, Vol. 79, No. 2, 2001, pp. 77{95 [32] Schaub, H., \Incorporating Secular Drifts into the Orbit Element Dierence Description of Relative Orbits," 13th AAS/AIAASpace Flight Mechanics Meeting, American Astronautical Society Paper 2003-115, 2003 [33] Sedwick, R. J., Miller, D. W., and Kong, E. M. C., "Mitigation of Dierential Perturbations in Clusters of Formation Flying Satellites," AAS/AIAA Space Flight Mechanics Meeting, American Astronautical Society, AAS Paper 99-124, Feb. 1999. [34] S. Schweighart and R. Sedwick, "High-Fidality Linearized J 2 Model for Satellite Formation Flight.", Journal of Guidance, Control and Dynamics, Vol.25, No.6, Nov-Dec 2002, pp1073- 1080 [35] T. Kasai, M. Oda and T. Suzuki, "Results of the ETS-7 Mission { Rendezvous Docking and Space Robotics Experiment", in 5th Int. Symp. on Articial Intelligence, Robotics and Auto. in Space, ESTEC/ESA, Nordwijk, The Netherlands, pp.299-306., 1999. [36] Tschauner, J. and Hempel, P., \Rendezvous Zu Einem In Elliptischer Bahn Umlaufenden Ziel," Acta Astronautica, Vol. 11, No. 2, 1965, pp. 104{109. [37] B. Udrea, and M. Nayak, "A Cooperative Multi-Satellite Mission for Controlled Active Debris Removal from Low Earth Orbit," in 2015 IEEE Aerospace Conference, 2015. 119 [38] Vaddi S.S., Vadali S. R., and Alfriend K.T., "Formation Flying: Accommodating Nonlinearity and Eccentricity Perturbations", Journal of Guidance Control and Dynamics, Vol. 26, No. 2, March-April 2013 [39] Vadali S.R., "Model for Linearized Satellite Relative Motion About a J 2 -Perturbed Mean Circular Orbit", Journal of Guidance Control and Dynamics, Vol. 32, No. 5, Sep.-Oct. 2009 [40] Vallado D.A., Fundamentals of Astrodynamics and Applications, 3rd Edition, Microcosm Press, 2007 120

Abstract (if available)

Abstract There has been an increasing interest in on-orbit autonomous servicing and repair of satellites as well as controlled active debris removal (ADR) in the space industry recently. One of the most challenging tasks for servicing/repair as well as for ADR is the rendezvous and docking with a non-cooperative tumbling resident space object (RSO). The work presented here deals with the optimal trajectory planning for proximity maneuvering and close approach for a nanosatellite with a tumbling RSO under high fidelity J₂ and quadratic drag perturbation model. ❧ The effect of J₂ and quadratic drag is studied on the relative motion and its importance as a requirement is established for implementing proximity operation algorithms for closed-loop attitude control and relative navigation. Current work implements a minimum-control-effort problem for a 3-D rendezvous to a tumbling object which considers a three-degree-of-freedom model that consists of six kinematics states, position, and velocity, for both the chaser/deputy and target/chief spacecraft. ❧ The control problem features additional path constraints with relative motion dynamics pertinent to the proximity space operations to match the perching state vector of the nanosatellite over a feature of interest of the RSO. The path to close approach, with a terminal constraint of a small but finite positive relative speed at contact is also discussed. The consequences of limited thrust and finite attitude maneuver time are taken into account and their effects on the closed-loop translation control of the nanosatellite are analyzed. Moreover, the research also elaborates on the homotopic nature of the optimal control trajectories for variations in different mission design parameters including inclination, maximum thrust, and altitude. ❧ A secondary control problem is formulated with LQR cost, this facilitates a full state feedback controller with user-defined state and control penalty. Optimal trajectories are studied under the defined penalty matrices. Finally, research shows a generic implementation of a neuromorphic controller that is powered by the population of neurons that has the ability to potentially do gain scheduling in real time.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Spacecraft trajectory optimization: multiple-impulse to time-optimal finite-burn conversion

PDF

Relative-motion trajectory generation and maintenance for multi-spacecraft swarms

PDF

Techniques for analysis and design of temporary capture and resonant motion in astrodynamics

PDF

Trajectory mission design and navigation for a space weather forecast

PDF

Application of the fundamental equation to celestial mechanics and astrodynamics

PDF

A declarative design approach to modeling traditional and non-traditional space systems

PDF

Optimization of relative motion rendezvous via modified particle swarm and primer vector theory

PDF

Optimizing element & system compliance of robotic, gecko adhesion-based grippers to unprepared space debris targets

PDF

New approaches to satellite formation-keeping and the inverse problem of the calculus of variations

PDF

Context-adaptive expandable-compact POMDPs for engineering complex systems

PDF

Optimization-based whole-body control and reactive planning for a torque controlled humanoid robot

Asset Metadata

Creator Patel, Parv (author)

Core Title Optimal guidance trajectories for proximity maneuvering and close approach with a tumbling resident space object under high fidelity J₂ and quadratic drag perturbation model

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Astronautical Engineering

Publication Date 11/15/2021

Defense Date 09/28/2021

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag astrodynamics,closed-loop,Drag,feedback,J₂,navigation,OAI-PMH Harvest,optimal control,orbital mechanics,perturbations

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Erwin, Daniel (committee chair), Gruntman, Mike (committee member), Madni, Azad (committee member), Udrea, Bogdan (committee member), Udwadia, Firdaus (committee member)

Creator Email parvpate@usc.edu,parvpatel@gmail.com

Permanent Link (DOI) https://doi.org/10.25549/usctheses-oUC17138475

Unique identifier UC17138475

Legacy Identifier etd-PatelParv-10221

Document Type Dissertation

Rights Patel, Parv

Type texts

Source 20211115-wayne-usctheses-batch-896-nissen (batch), University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the author, as the original true and official version of the work, but does not grant the reader permission to use the work if the desired use is covered by copyright. It is the author, as rights holder, who must provide use permission if such use is covered by copyright. The original signature page accompanying the original submission of the work to the USC Libraries is retained by the USC Libraries and a copy of it may be obtained by authorized requesters contacting the repository e-mail address given.

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA

Repository Email cisadmin@lib.usc.edu