Attitude Trajectory Optimization for Agile Satellites in Autonomous Remote Sensing Constellation
Emmanuel Sin, Sreeja Nag, Vinay Ravindra, Alan Li, Murat Arcak
AAttitude Trajectory Optimization for Agile Satellites inAutonomous Remote Sensing Constellations
Emmanuel Sin ∗ , Murat Arcak † University of California, Berkeley, California, 94720, U.S.A.
Sreeja Nag ‡ , Vinay Ravindra § , Alan Li ¶ NASA Ames Research Center, Moffett Field, California, 94035, U.S.A.
Agile attitude maneuvering maximizes the utility of remote sensing satellite constellations.By taking into account a satellite’s physical properties and its actuator specifications, we mayleverage the full performance potential of the attitude control system to conduct agile remotesensing beyond conventional slew-and-stabilize maneuvers. Employing a constellation of agilesatellites, coordinated by an autonomous and responsive scheduler, can significantly increaseoverall response rate, revisit time and global coverage for the mission. In this paper, weuse recent advances in sequential convex programming (SCP) based trajectory optimizationto enable rapid-target acquisition, pointing and tracking capabilities for a scheduler-basedconstellation. We present two problem formulations. The
Minimum-Time Slew Optimal ControlProblem determines the minimum time, required energy, and optimal trajectory to slew betweenany two orientations given nonlinear quaternion kinematics, gyrostat and actuator dynamics,and state/input constraints. By gridding the space of 3D rotations and efficiently solving thisproblem on the grid, we produce lookup tables or parametric fits off-line that can then be usedon-line by a scheduler to compute accurate estimates of minimum-time and maneuver energyfor real-time constellation scheduling. The estimates allow an optimization-based scheduler toproduce target-remote-sensing and data-downlinking schedules that are dynamically feasiblefor each satellite and optimal for the constellation. The
Minimum-Effort Multi-Target PointingOptimal Control Problem is used on-line by each satellite to produce continuous attitude-stateand control-input trajectories that realize a given schedule while minimizing attitude errorand control effort. The optimal trajectory may then be achieved by a low-level trackingcontroller. This onboard trajectory generation and tracking scheme is possible due to real-time, efficient SCP implementations. We demonstrate our approach with a numerical examplethat uses simulation data for a reference satellite in Sun-synchronous orbit passing over globally-distributed, Earth-observation targets.
I. Nomenclature 𝐽 = mass moment of inertia matrix of rigid body in body-fixed frame, [kg · m ] 𝐽 𝑟 = mass moment of inertia of actuation rotor about its axis of rotation, [kg · m ] 𝐴 𝑟 = actuator Jacobian, each column representing a rotor’s axis of rotation w.r.t. body-fixed frame q = unit quaternion parameterizing attitude of rigid body’s body-fixed frame w.r.t. inertial frame, [ ] 𝝎 = angular velocity of rigid body about body-fixed frame axes, [rad/s] r i = angular momentum of rotor 𝑖 about its axis, [Nms] u 𝑖 = torque of rotor 𝑖 about its rotor axis, [N] 𝑡 = time, [s] 𝜏 = normalized time, [ ] ∗ Graduate Student, Department of Mechanical Engineering, [email protected], AIAA Student Member. † Professor, Department of Electrical Engineering and Computer Sciences, [email protected]. ‡ Senior Research Scientist, Bay Area Environmental Research Institute, [email protected], AIAA Senior Member. § Research Scientist, Bay Area Environmental Research Institute, [email protected], AIAA Member. ¶ Research Scientist, Bay Area Environmental Research Institute, [email protected]. a r X i v : . [ ee ss . S Y ] F e b I. Introduction
Due to the proliferation of launch providers to low Earth orbit (LEO) and the trend towards smaller, cost-efficientspacecraft, satellite constellations are enabling scientific missions and commercial applications that are otherwiseimpossible with a single, larger satellite. For example, a constellation of satellites in LEO may be coordinated to pointtowards coastal regions around the world to measure ocean color, atmospheric properties, phytoplankton concentrations,and ultimately assess the health of global coral reef ecosystems [1]. LEO constellations may also be employed tomeasure episodic precipitation and subsequent water flow in flood-prone cities [2]-[3]. Furthermore, constellations maybe tasked to measure soil moisture in targeted regions to assess the risk of wildfires [4]. In addition to climate andenvironment monitoring, Earth-observing constellations are providing commercial and economic value by measuring,for instance, agricultural crop growth, infrastructure development, or logistical activity at airports and shipping containerroutes [5]. Beyond these Earth observation applications, LEO satellite constellations are also enabling space-basedInternet and telecommunication services [6].To enhance performance as an interconnected system, the satellites of a constellation can be precisely coordinatedthroughout each stage of a mission. Once a group of satellites are deployed into a desired orbital plane by a launchvehicle, the satellites enter the orbit acquisition stage where they must be phased relative to each other to achieve desiredangular spacing. In [7], it is assumed that the orientation of the satellites, with respect to their orbital velocities, may becontrolled so that either a minimum or maximum surface area is exposed to the atmosphere that exists in low Earthorbit. By inducing either a low or high atmospheric drag on each satellite using a bang-bang control approach, theresulting differential drag between satellites is used to separate them into an equally-spaced formation. In [8], simulatedannealing is used to design time-optimal differential drag commands for a group of up to 100 satellites and the method isdemonstrated on an actual constellation deployed in LEO. Assuming accurate attitude pointing, [9] formulates a linearprogram that produces differential drag commands taking on continuous values between the minimum and maximumvalues, allowing the constellation to form not only under a minimum-time objective but also with a maximum altitude (orequivalently, maximum constellation lifetime) objective. A distributed controller approach is presented in [10], where itis assumed that the attitude of the satellites can be controlled to apply continuous low-thrust in the appropriate direction.Regardless of where the control authority is derived from (e.g., differential drag commands or thrust commands), theseconstellation acquisition methods require accurate attitude pointing.In the subsequent operational stage , the satellites must perform various scheduled tasks including targeted remotesensing (e.g., imaging, radiometry), downlinking/uplinking data to/from ground stations, orbital station-keeping andother maintenance activities. In [1]-[4], an automated scheduler is developed to run autonomously on either groundstations (with schedules uplinked to satellites) or onboard in a distributed approach. Based on dynamic programming ormixed integer programming, the scheduler produces imaging schedules for each satellite that maximizes the numberof observations and/or observation times for the constellation as a whole. In addition to maximizing spatial imagingcoverage, [5] also addresses the problem of maximizing data downlinked by the constellation. We note that theschedulers in these works make inherent assumptions on the agility and pointing performance of the attitude controlsubsystem on board each satellite. For instance, in order to produce a feasible schedule that provides sufficient time toslew between desired orientations, the scheduler must know the dynamically-feasible, minimum slew time betweenany two desired orientations. This minimum slew time depends on the physical properties of the spacecraft (i.e., massmoment of inertia, actuator configuration) and actuator constraints (e.g., maximum rotor torque and momentum). Oncean optimized schedule is produced for the constellation, each satellite must then generate and track an attitude trajectorythat realizes its given schedule. Such an attitude trajectory may be optimized to minimize desired-actual attitude error atspecific times of the trajectory and minimize control effort over the course of the trajectory.Predicting that future spacecraft will require agile attitude control systems that provide rapid multi-target pointingand tracking capabilities, [11] proposes a feedback regulator to conduct large-angle, rest-to-rest slew maneuvers usingthe Euler’s eigenaxis rotation between any two orientations. Resembling feedback linearization, the proposed lawintroduces a nonlinear term to cancel out the coupling between body angular velocities and replaces it with linearerror-quaternion and body rate feedback terms. In [12], a large class of attitude tracking control laws that have thegeneral form of proportional-derivative (PD) feedback and feedforward compensation is obtained with proofs of theirglobal asymptotic stability in the closed-loop. Minimum-time, rest-to-rest slew maneuvers for an inertially symmetricrigid spacecraft with independent three-axis controls are studied in [13], which shows that the optimal maneuver isnot, in general, an eigenaxis rotation but one that includes significant nutation of the instantaneous axis of rotation.Furthermore, the structure of the optimal control is different for small and large reorientation angles. While [11]-[13]focus on rigid body dynamics under ideal, body-fixed control torques, [14] considers actuator dynamics and presents afeedback control logic that produces near minimum-time eigenaxis slew maneuvers under actuator saturation, slew2ate limit, and control bandwidth limit. The control laws are used in [15] to demonstrate rapid multi-ground-targetacquisition by stepping through reference set-points that define successive scan trajectories. We note that the attitudecontrol strategy in [1] also uses a minimum-time, eigenaxis slew maneuver and switches to a PD control law for smallangles. Based on closed-loop simulations with this control strategy, a polynomial fit of minimum maneuver time as afunction of the eigenaxis slew angle is used in a scheduler [2].Although eigenaxis-based minimum-time control laws may be applied to general minimum-time problems to producenear-optimal maneuvers, they do not explicitly address nonlinearities that arise from, for instance, actuator dynamicsnor do they explicitly consider general state and input constraints (e.g., bounds on momentum, power and energy).Moreover, when directly applying a feedback control strategy to track a desired sequence of discrete orientations, thefeedback gains between each pair of orientations must be carefully tuned to achieve settling times that completely satisfya desired pointing schedule, without missing any targets. In these minimum-time or minimum-attitude-error situations,we may use trajectory optimization methods to not only explicitly deal with nonlinearities and constraints, but to alsoautomate the generation of continuous attitude trajectories for autonomous execution.As surveyed by Betts [16], there is an expanse of literature on trajectory optimization methods that may becharacterized by, for example, the solution approach (i.e., indirectly satisfying necessary conditions for optimality ordirectly solving a transcribed version of the problem), or the transcription process (shooting versus collocation). In [17],various direct collocation methods are introduced while [18] focuses on an indirect pseudospectral method that has beenused in practice, for example, to control the International Space Station with a zero-propellant maneuver. Regardless ofthe approach, many trajectory optimization methods treat a problem in its original nonlinear, non-convex form, requiringthe use of a nonlinear programming solver [19], which may be computationally inefficient with long solution times.Recent advances in sequential convex programming (SCP) have enabled efficient computation of locally optimaltrajectories for nonlinear systems with non-convex constraints and objectives. SCP is an iterative method that repeatedlyformulates and solves a convex, finite-dimensional parameter optimization problem that approximates the originalnon-convex optimal control problem. A convex formulation is typically achieved by linearizing the nonlinear systemaround a nominal trajectory (i.e., the solution from the previous iteration) and approximating any non-convex constraintsand objective with Taylor series expansions. Fast and reliable Interior Point Method algorithms [20] may be used to solvethese convex subproblems. In the early works of [21] and [22], successive convexification, a specific implementation ofthe SCP method, is introduced to find minimum-fuel and minimum-time trajectories in a 6-DOF rocket landing problem.In [23]-[26], certain details of the implementation, including choice of discretization method, constraint satisfactionbetween temporal nodes, convexification of non-convex constraints, and algorithm convergence properties are explored.Finally, the works of [27]-[33] introduce state-triggered constraints and address real-time, onboard implementations thatmay produce solutions in a fraction of a second.
A. Main Contributions
The main contributions of this work are the formulation and application of two optimal control problems (OCPs):1) An off-line method to produce accurate estimates of the optimal slew time and required energy for a referencesatellite to maneuver between any two arbitrary orientations. The estimates may be used by a scheduler toproduce dynamically-feasible pointing schedules for each satellite of a constellation. We call the associatedproblem:
Minimum-Time Slew OCP .2) An on-line method to produce optimal attitude trajectories that satisfy desired multi-target pointing scheduleswith minimum control effort. This method may be run on-board to generate trajectories for a low-level trackingcontroller to then follow. The corresponding problem is called:
Minimum-Effort Multi-Target Pointing OCP .These two contributions can be applied together to support the scheduler. The first contribution allows a constellationscheduler to make informed multi-target pointing decisions (i.e., pointing schedules) based on an accurate assessmentof a spacecraft’s maneuvering capabilities. The second contribution allows the spacecraft to plan feasible attitudetrajectories that satisfy the desired pointing schedules generated by the scheduler.
B. Organization of Paper
After introducing the dynamical model and attitude parameterization conventions assumed in this paper, we presenttwo problem formulations with minimum-time slew and minimum-effort multi-target pointing objectives. We describethe process of transcribing the continuous-time trajectory optimization problems into convex, finite-dimensionalparameter optimization problems. We then review state-of-the-art techniques in sequential convex programming (SCP)that are used to approximately solve the original problems. Finally, we apply the formulations to numerical examples.3 . Preliminaries
To study attitude motion, a satellite may be modeled as a gyrostat [34], consisting of the platform (i.e., spacecraftbus) and actuation rotors (e.g., momentum/reaction wheels). The state vector for this system of rigid bodies may berepresented as: x = [ q (cid:62) 𝝎 (cid:62) r (cid:62) ] (cid:62) = [ 𝑞 𝑞 𝑞 𝑞 𝜔 𝑥 𝜔 𝑦 𝜔 𝑧 𝑟 𝑟 𝑟 𝑟 ] (cid:62) (1)where q is the unit quaternion (i.e., (cid:107) q (cid:107) =
1) that describes the orientation of the spacecraft’s body fixed frame withrespect to an inertial frame, 𝝎 is the angular velocity vector of the spacecraft with components expressed about thebody-fixed axes, and r represents the angular momentum of each spinning rotor about its axis of rotation (we assume thespacecraft has four actuation rotors). The input vector consists of the torques produced by each rotor about its axis: u = [ 𝑢 𝑢 𝑢 𝑢 ] (cid:62) (2)The differential equations describing the motion of this 11-state, 4-input system consist of the quaternion kinematics,gyrostat equation, and a single-integrator model for the rotors: (cid:164) x : = 𝑑𝑑𝑡 x = f ( x , u ) : = Ω q − 𝐽 − (cid:0) 𝝎 × ( 𝐽 𝝎 + 𝐴 𝑟 r ) + 𝐴 𝑟 u (cid:1) u (3)The positive definite matrix 𝐽 represents the mass moment of inertia of the spacecraft in the body-fixed frame, 𝐴 𝑟 ∈ R × is the actuator Jacobian where each column represents a rotor’s axis of rotation with respect to the body-fixed frame,and Ω and 𝝎 × are skew-symmetric matrices defined as: Ω : = 𝜔 𝑧 - 𝜔 𝑦 𝜔 𝑥 - 𝜔 𝑧 𝜔 𝑥 𝜔 𝑦 𝜔 𝑦 - 𝜔 𝑥 𝜔 𝑧 - 𝜔 𝑥 - 𝜔 𝑦 - 𝜔 𝑧 𝝎 × : = 𝜔 𝑧 𝜔 𝑦 𝜔 𝑧 𝜔 𝑥 - 𝜔 𝑦 𝜔 𝑥 (4)Note that we use the notation [ ] × to map a vector into a corresponding skew symmetric matrix.In this paper, we follow the unit quaternion (or Euler Symmetric Parameters) convention used in [35] and [36] forattitude parameterization, where the vector part is stacked on top of the scalar part: q : = [ q (cid:62) 𝑣 𝑞 𝑠 ] (cid:62) : = [ 𝑞 𝑞 𝑞 𝑞 ] (cid:62) (5)Furthermore, we denote the quaternion conjugate and the identity quaternion, respectively, as: q + : = [ - q (cid:62) 𝑣 𝑞 𝑠 ] (cid:62) = [ - 𝑞 - 𝑞 - 𝑞 𝑞 ] (cid:62) (6) q I : = [ (cid:62) ] (cid:62) = [ ] (cid:62) (7)As defined in [36] and used in [11] and [14], we may measure the attitude error between a given quaternion q and adesired quaternion ¯ q by computing the error-quaternion: q 𝑒 = 𝑞 𝑒 𝑞 𝑒 𝑞 𝑒 𝑞 𝑒 : = ¯ q + q = ¯ 𝑞 ¯ 𝑞 - ¯ 𝑞 - ¯ 𝑞 - ¯ 𝑞 ¯ 𝑞 ¯ 𝑞 - ¯ 𝑞 ¯ 𝑞 - ¯ 𝑞 ¯ 𝑞 - ¯ 𝑞 ¯ 𝑞 ¯ 𝑞 ¯ 𝑞 ¯ 𝑞 𝑞 𝑞 𝑞 𝑞 (8)We note that the Hamilton product (or any quaternion multiplication operation) between two quaternions is associativeand distributive, but not commutative in general, i.e., ¯ q + q ≠ q ¯ q + . III. Problem Formulations
In this section we introduce two optimal control problem (OCP) formulations to be solved with sequential convexprogramming and applied to numerical examples involving a rotor-actuated satellite under agile-pointing operations. Asstated, these two formulations have the same constraints and differ only in the objective and decision variables. Weconclude this section by discussing how the objectives and constraints may be modified depending on the application.4 . Minimum-Time Slew
The
Minimum-Time Slew OCP addresses the time-optimal, large-angle, rest-to-rest slew problem. The integralterm in (9) captures the minimum-time objective and we note that the decision variables consist of the gyrostat’s stateand input as well as the final time 𝑡 𝑓 . The gyrostat equation (11) may be modified to include other actuators, such asmagnetorquers and propulsive thrusters. We may also include modeled environmental disturbances, including momentsdue to gravity gradient, atmospheric drag, solar radiation pressure, and the magnetic field of Earth. Furthermore, theactuator dynamics in (12) may use higher-fidelity models that consider, for example, brushless DC motor dynamicsand rotational friction. Maximum rotor momentum and maximum rotor torque bounds are enforced by (13) - (14),respectively. Finally, the initial and final conditions are given in (15) - (16), where the overbar notation signifies adesired or given quantity.minimize q , 𝝎 , r , u ,𝑡 𝑓 ∫ 𝑡 𝑓 𝑡 𝑑𝑡 (9)s.t. (cid:164) q ( 𝑡 ) = Ω ( 𝑡 ) q ( 𝑡 ) ∀ 𝑡 ∈ [ 𝑡 , 𝑡 𝑓 ] (10) (cid:164) 𝝎 ( 𝑡 ) = − 𝐽 − (cid:0) 𝝎 × ( 𝑡 ) ( 𝐽 𝝎 ( 𝑡 ) + 𝐴 𝑟 r ( 𝑡 )) + 𝐴 𝑟 u ( 𝑡 ) (cid:1) ∀ 𝑡 ∈ [ 𝑡 , 𝑡 𝑓 ] (11) (cid:164) r ( 𝑡 ) = u ( 𝑡 ) ∀ 𝑡 ∈ [ 𝑡 , 𝑡 𝑓 ] (12) 𝑟 𝑖 ( 𝑡 ) < 𝑟 𝑚𝑎𝑥 ∀ 𝑖 = , . . . , ∀ 𝑡 ∈ [ 𝑡 , 𝑡 𝑓 ] (13) 𝑢 𝑖 ( 𝑡 ) < 𝑢 𝑚𝑎𝑥 ∀ 𝑖 = , . . . , ∀ 𝑡 ∈ [ 𝑡 , 𝑡 𝑓 ] (14) q ( 𝑡 ) = ¯ q 𝑖𝑛𝑖𝑡𝑎𝑙 , 𝝎 ( 𝑡 ) = (15) q ( 𝑡 𝑓 ) = ¯ q 𝑓 𝑖𝑛𝑎𝑙 , 𝝎 ( 𝑡 𝑓 ) = (16) B. Minimum-Effort Multi-Target Pointing
In the following
Minimum-Effort Multi-Target Pointing OCP , we maintain the same constraints as the problem abovebut change the objective. Furthermore, the final time 𝑡 𝑓 is no longer a decision variable but a fixed parameter. Given aschedule (i.e., sequence) of desired quaternions and angular velocities at specified times, { 𝑡 𝑘 , ¯ q ( 𝑡 𝑘 ) , ¯ 𝝎 ( 𝑡 𝑘 )} ∀ 𝑘 ∈ K ,where K is a finite set of indices, the objective in (17) minimizes the error at those discrete time points while alsominimizing the continuous control effort. The control effort used during the trajectory may be expressed as the energyof input signal, defined as ∫ 𝑡 𝑓 𝑡 (cid:107) u ( 𝑡 )(cid:107) 𝑑𝑡 . We note that minimizing the last term in 17 is equivalent to minimizing theenergy of the input signal. The user-defined parameter 𝛾 > 𝜌 > q , 𝝎 , r , u ∑︁ 𝑘 ∈ K (cid:8)(cid:13)(cid:13) ¯ q + ( 𝑡 𝑘 ) q ( 𝑡 𝑘 ) − q I (cid:13)(cid:13) + 𝛾 (cid:107) 𝝎 ( 𝑡 𝑘 ) − ¯ 𝝎 ( 𝑡 𝑘 )(cid:107) (cid:9) + 𝜌 ∫ 𝑡 𝑓 𝑡 (cid:107) u ( 𝑡 )(cid:107) 𝑑𝑡 (17)s.t. (cid:164) q ( 𝑡 ) = Ω ( 𝑡 ) q ( 𝑡 ) ∀ 𝑡 ∈ [ 𝑡 , 𝑡 𝑓 ] (18) (cid:164) 𝝎 ( 𝑡 ) = − 𝐽 − (cid:0) 𝝎 × ( 𝑡 ) ( 𝐽 𝝎 ( 𝑡 ) + 𝐴 𝑟 r ( 𝑡 )) + 𝐴 𝑟 u ( 𝑡 ) (cid:1) ∀ 𝑡 ∈ [ 𝑡 , 𝑡 𝑓 ] (19) (cid:164) r ( 𝑡 ) = u ( 𝑡 ) ∀ 𝑡 ∈ [ 𝑡 , 𝑡 𝑓 ] (20) 𝑟 𝑖 ( 𝑡 ) < 𝑟 𝑚𝑎𝑥 ∀ 𝑖 = , . . . , ∀ 𝑡 ∈ [ 𝑡 , 𝑡 𝑓 ] (21) 𝑢 𝑖 ( 𝑡 ) < 𝑢 𝑚𝑎𝑥 ∀ 𝑖 = , . . . , ∀ 𝑡 ∈ [ 𝑡 , 𝑡 𝑓 ] (22) q ( 𝑡 ) = ¯ q 𝑖𝑛𝑖𝑡𝑎𝑙 , 𝝎 ( 𝑡 ) = (23) q ( 𝑡 𝑓 ) = ¯ q 𝑓 𝑖𝑛𝑎𝑙 , 𝝎 ( 𝑡 𝑓 ) = (24)We note that an alternative “constraint formulation” of the Minimum-Effort Multi-Target Pointing OCP is to remove theattitude penalty terms in the objective and implement them as constraints: (cid:13)(cid:13) ¯ q + ( 𝑡 𝑘 ) q ( 𝑡 𝑘 ) − q I (cid:13)(cid:13) ≤ 𝜖 q ∀ 𝑘 ∈ K (25) (cid:107) 𝝎 ( 𝑡 𝑘 ) − ¯ 𝝎 ( 𝑡 𝑘 )(cid:107) ≤ 𝜖 𝝎 ∀ 𝑘 ∈ K (26)5here 𝜖 q ≥ 𝜖 𝝎 ≥ 𝑃 𝑚𝑎𝑥 and maximum energy consumed 𝐸 𝑚𝑎𝑥 by the attitude control system (i.e., all four rotors combined): ∑︁ 𝑖 = (cid:12)(cid:12)(cid:12)(cid:12) 𝑢 𝑖 ( 𝑡 ) · 𝐽 𝑟 𝑟 𝑖 ( 𝑡 ) (cid:12)(cid:12)(cid:12)(cid:12) < 𝑃 𝑚𝑎𝑥 ∀ 𝑡 ∈ [ 𝑡 , 𝑡 𝑓 ] (27) ∫ 𝑡 𝑓 𝑡 (cid:40) ∑︁ 𝑖 = (cid:12)(cid:12)(cid:12)(cid:12) 𝑢 𝑖 ( 𝑡 ) · 𝐽 𝑟 𝑟 𝑖 ( 𝑡 ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:41) 𝑑𝑡 < 𝐸 𝑚𝑎𝑥 (28)where 𝐽 𝑟 is the rotor inertia. Apart from the equations of motion, we note that the instantaneous power and energyconstraints are non-convex due to being bilinear in the decision variables of u and r . In the following section we discusshow such non-convex constraints can be approximated as convex constraints. IV. Trajectory Optimization
In this section, we first describe how a general optimal control problem (OCP) is transcribed into a convex,finite-dimensional parameter optimization problem (OPT). We then review the sequential convex programming methodused in this paper. Consider the following continuous-time dynamical system: (cid:164) x ( 𝑡 ) : = 𝑑𝑑𝑡 x ( 𝑡 ) = f ( x ( 𝑡 ) , u ( 𝑡 )) ∀ 𝑡 ∈ [ 𝑡 , 𝑡 𝑓 ] (29)defined over the given time span, where x ( 𝑡 ) ∈ R 𝑛 𝑥 is the state of the system and u ( 𝑡 ) ∈ R 𝑛 𝑢 is the input to the system. A. Time Normalization
In [22], a procedure is introduced to transform a free-final-time optimal control problem into a finite-dimensionalparameter optimization problem by normalizing time in the dynamical system (29) by a normalization factor 𝑡 𝑓 : 𝜏 : = 𝑡𝑡 𝑓 (30)We treat 𝑡 𝑓 as a decision variable in a free-final-time problem (e.g., minimum-time OCP) and as a constant parameterin a fixed-final-time problem (e.g., minimum-attitude-error-and-control-effort OCP). We now express the time spanof the dynamical system in terms of the normalized time: 𝑡 𝑡 𝑓 ≤ 𝜏 ≤
1. Since 𝑡 = 𝑡 𝑓 𝜏 ⇒ 𝑑𝑡 = 𝑡 𝑓 𝑑𝜏 ⇒ 𝑑𝑡𝑑𝜏 = 𝑡 𝑓 , thederivative of the scaled state with respect to normalized time is: x (cid:48) ( 𝑡 ) : = 𝑑𝑑𝜏 x ( 𝑡 ) = 𝑑𝑡𝑑𝜏 𝑑𝑑𝑡 x ( 𝑡 ) = 𝑡 𝑓 f ( x ( 𝑡 ) , u ( 𝑡 )) ∀ 𝑡 ∈ [ 𝑡 , 𝑡 𝑓 ] (31)Since 𝑡 can be expressed as a function of 𝜏 and assuming 𝑡 =
0, we can represent (31) in terms of normalized time: x (cid:48) ( 𝑡 ( 𝜏 )) = x (cid:48) ( 𝜏 ) = 𝑡 𝑓 f ( x ( 𝜏 ) , u ( 𝜏 )) = : F ( x ( 𝜏 ) , u ( 𝜏 ) , 𝑡 𝑓 ) ∀ 𝑡 ∈ [ , ] (32)In this final expression, it is clear that the dynamical system is a function of state x , input u , and final time 𝑡 𝑓 . B. Linearization
Assuming the dynamical system described by (32) is non-convex but differentiable, it can be approximated as convexin x , u , and 𝑡 𝑓 with a first-order Taylor expansion about a given trajectory: x (cid:48) ( 𝜏 ) = F ( x ( 𝜏 ) , u ( 𝜏 ) , 𝑡 𝑓 ) ≈ 𝐴 ( 𝜏 ) x ( 𝜏 ) + 𝐵 ( 𝜏 ) u ( 𝜏 ) + Σ ( 𝜏 ) 𝑡 𝑓 + e ( 𝜏 ) ∀ 𝑡 ∈ [ , ] (33)where we denote the first-order partial derivative matrices of F ( x ( 𝜏 ) , u ( 𝜏 ) , 𝑡 𝑓 ) evaluated about { ¯ x ( 𝜏 ) , ¯ u ( 𝜏 ) , ¯ 𝑡 𝑓 } as 𝐴 ( 𝜏 ) : = 𝐷 x F ( ¯ x ( 𝜏 ) , ¯ u ( 𝜏 ) , ¯ 𝑡 𝑓 ) (34) 𝐵 ( 𝜏 ) : = 𝐷 u F ( ¯ x ( 𝜏 ) , ¯ u ( 𝜏 ) , ¯ 𝑡 𝑓 ) (35) Σ ( 𝜏 ) : = 𝐷 𝑡 𝑓 F ( ¯ x ( 𝜏 ) , ¯ u ( 𝜏 ) , ¯ 𝑡 𝑓 ) = f ( ¯ x ( 𝜏 ) , ¯ u ( 𝜏 )) (36)6ith the following dynamical approximation offset term: e ( 𝜏 ) : = −( 𝐴 ( 𝜏 ) ¯ x ( 𝜏 ) + 𝐵 ( 𝜏 ) ¯ u ( 𝜏 )) (37) C. Discretization
As described in [37], the process of converting an optimal control problem into a parameter optimization problembegins by dividing the time duration of the optimal control problem into intervals using 𝐾 temporal nodes. The nodesmay be chosen to be equally spaced, creating 𝐾 − = : 𝜏 < 𝜏 < · · · < 𝜏 𝑘 < . . . < 𝜏 𝐾 − < 𝜏 𝐾 : = x 𝑘 : = x ( 𝜏 𝑘 ) , u 𝑘 : = u ( 𝜏 𝑘 ) , e 𝑘 : = e ( 𝜏 𝑘 ) ∀ 𝑘 = , . . . , 𝐾 (39)For use in the following section, we collectively refer to the decision variables that we have influence over as: z 𝑘 : = [ x (cid:62) 𝑘 , u (cid:62) 𝑘 , 𝑡 𝑓 ] (cid:62) ∀ 𝑘 = , . . . , 𝐾 (40)We use the First-Order Hold (FOH)-interpolation based discretization method in our implementation. As demonstratedin [24], the FOH discretization provides fast computational time and achieves similar accuracy when compared to moreadvanced pseudospectral methods. Furthermore, it was shown that if convex input constraints are satisfied at the nodes,then inter-nodal convex input constraint satisfaction is also guaranteed [23]. For convenience, we review the processused in [21]-[32] below.The FOH interpolation represents the input within each of the 𝐾 − u ( 𝜏 ) = 𝜆 − 𝑘 u 𝑘 + 𝜆 + 𝑘 u 𝑘 + ∀ 𝜏 ∈ [ 𝜏 𝑘 , 𝜏 𝑘 + ] , 𝑘 = , . . . , 𝐾 − 𝜆 − 𝑘 : = 𝜏 𝑘 + − 𝜏𝜏 𝑘 + − 𝜏 𝑘 , 𝜆 + 𝑘 : = 𝜏 − 𝜏 𝑘 𝜏 𝑘 + − 𝜏 𝑘 (42)The exact discretization of (33) is then: x 𝑘 + = 𝐴 𝑘 x 𝑘 + 𝐵 − 𝑘 u 𝑘 + 𝐵 + 𝑘 u 𝑘 + + Σ 𝑘 𝑡 𝑓 + e 𝑘 ∀ 𝑘 = , . . . , 𝐾 − 𝐴 𝑘 : = Φ ( 𝜏 𝑘 + , 𝜏 𝑘 ) (44) 𝐵 − 𝑘 : = 𝐴 𝑘 ∫ 𝜏 𝑘 + 𝜏 𝑘 Φ − ( 𝜏, 𝜏 𝑘 ) 𝜆 − 𝑘 ( 𝜏 ) 𝐵 ( 𝜏 ) 𝑑𝜏 (45) 𝐵 + 𝑘 : = 𝐴 𝑘 ∫ 𝜏 𝑘 + 𝜏 𝑘 Φ − ( 𝜏, 𝜏 𝑘 ) 𝜆 + 𝑘 ( 𝜏 ) 𝐵 ( 𝜏 ) 𝑑𝜏 (46) Σ 𝑘 : = 𝐴 𝑘 ∫ 𝜏 𝑘 + 𝜏 𝑘 Φ − ( 𝜏, 𝜏 𝑘 ) Σ ( 𝜏 ) 𝑑𝜏 (47) e 𝑘 : = 𝐴 𝑘 ∫ 𝜏 𝑘 + 𝜏 𝑘 Φ − ( 𝜏, 𝜏 𝑘 ) w ( 𝜏 ) 𝑑𝜏 (48)The state transition matrix satisfies the following differential equation and initial condition within each interval: 𝑑𝑑𝜏 Φ ( 𝜏, 𝜏 𝑘 ) = 𝐴 ( 𝜏 ) Φ ( 𝜏, 𝜏 𝑘 ) , Φ ( 𝜏 𝑘 , 𝜏 𝑘 ) = 𝐼 𝑛 𝑥 (49)In practice, the integrands of (45)-(48) along with (49) and (32) are numerically integrated from the start to the end ofeach interval using a nominal trajectory { ¯ x 𝑘 , ¯ u 𝑘 , ¯ 𝑡 𝑓 , ¯ e 𝑘 } for 𝑘 = , . . . , 𝐾 −
1. Note that since we initialize the numericalintegration for each interval with points from a nominal trajectory, rather than the terminal points found by integrationin the previous temporal interval, this approach resembles a multiple-shooting discretization method, which has shownto improve convergence of the SCP algorithm by keeping solutions closer to the nominal trajectory. In contrast, a singleshooting method allows approximation errors to grow in later temporal intervals [24].7 . Sequential Convex Programming Method
For the example in this paper, we use the Penalized Trust Region (PTR) variant of SCP described in [33]. A keydifference with Successive Convexification (SCvx) studied in [26] is that PTR treats trust regions as soft constraintsplaced in the objective whereas SCvx enforces hard trust region constraints that are updated based on a rule. Anadvantage of SCvx is that convergence of this method is guaranteed. However, the method employs slack variables thatmay cause the approximately solved problem to be far from the original problem if they take on non-zero values inthe solution. Hence, a converged solution may be infeasible for the original problem. In the following subsections wedescribe the PTR implementation of virtual controls and trust regions, and an approach to constraint convexification.
1. Virtual Controls
While executing sequential convex programming on a trajectory optimization problem, the approximated convexproblem may become infeasible. This artificial infeasibility [22] is frequently encountered in the early iterations of thealgorithm when the dynamics are linearized about a poor initial guess. To alleviate this issue, slack variables called virtual controls are added to the discrete-time equations of motion: (43) x 𝑘 + = 𝐴 𝑘 x 𝑘 + 𝐵 − 𝑘 u 𝑘 + 𝐵 + 𝑘 u 𝑘 + + Σ 𝑘 𝑡 𝑓 + e 𝑘 + v 𝑘 ∀ 𝑘 = , . . . , 𝐾 − 𝐽 𝑣𝑐 = 𝑤 𝑣𝑐 𝐾 ∑︁ 𝑘 = (cid:107) v 𝑘 (cid:107) (51)where 𝑤 𝑣𝑐 is a large positive weight. Minimization of the 1-norm term encourages sparsity in the virtual control vector.
2. Trust Regions
To ensure that the solver does not stray too far from a nominal trajectory where the linearized model becomes lessaccurate, we implement a trust region cost term. In PTR, the deviation of decision variables from the solution of theprevious iteration is penalized with the 2-norm of weighted deviations: 𝐽 𝑡𝑟 = 𝐾 ∑︁ 𝑘 = (cid:107) 𝑊 𝑡𝑟 ( z 𝑘 − ¯ z 𝑘 )(cid:107) (52)Assuming the system has been scaled, the weight matrix 𝑊 𝑡𝑟 may be designed as a diagonal matrix with 𝑤 𝑡𝑟 > 𝑊 𝑡𝑟 : = 𝑤 𝑡𝑟 · 𝐼 𝐾 ( 𝑛 𝑥 + 𝑛 𝑢 + ) (53)
3. Constraint Convexification
Let us consider a general (non-convex) constraint on the continuous-time variables: g ( x ( 𝜏 ) , u ( 𝜏 ) , 𝑡 𝑓 ) ≤ ∀ 𝜏 ∈ [ , ] (54)This constraint can be approximated by linearizing it about a nominal trajectory:˜ g ( x ( 𝜏 ) , u ( 𝜏 ) , 𝑡 𝑓 ) : = g ( ¯ x ( 𝜏 ) , ¯ u ( 𝜏 ) , ¯ 𝑡 𝑓 ) + 𝐷 g ( ¯ x ( 𝜏 ) , ¯ u ( 𝜏 ) , ¯ 𝑡 𝑓 ) (cid:169)(cid:173)(cid:173)(cid:171) x ( 𝜏 ) u ( 𝜏 ) 𝑡 𝑓 − ¯ x ( 𝜏 ) ¯ u ( 𝜏 ) ¯ 𝑡 𝑓 (cid:170)(cid:174)(cid:174)(cid:172) ≤ ∀ 𝜏 ∈ [ , ] (55)We further approximate the constraint by explicitly enforcing it only at the time nodes of our discretization:˜ g ( x 𝑘 , u 𝑘 , 𝑡 𝑓 ) = ˜ g ( z 𝑘 ) ≤ ∀ 𝑘 = , . . . , 𝐾 (56)Inter-nodal constraint satisfaction is not guaranteed. However, the constraints at the nodes can be carefully designed sothat constraints are enforced for all time [23]. In a similar fashion, any non-convex cost terms may also be approximatedwith a Taylor series expansion about a nominal trajectory. 8e may also approximate constraints involving integral terms. For example: ∫ 𝑡 𝑓 𝑡 ℎ ( x ( 𝑡 ) , u ( 𝑡 )) 𝑑𝑡 ≤ = ⇒ 𝑡 𝑓 ∫ ℎ ( x ( 𝜏 ) , u ( 𝜏 )) 𝑑𝜏 ≤ = ⇒ 𝐾 − 𝐾 ∑︁ 𝑘 = 𝑡 𝑓 ℎ ( x 𝑘 , u 𝑘 ) ≤
4. Transcribed Convex Program
The
Minimum-Time OCP , described by equations (9) - (16), may be transcribed into the following convex,finite-dimensional parameter optimization problem that we call
Minimum-Time Slew OPT :minimize { z 𝑘 , v 𝑘 } 𝐾𝑘 = 𝑡 𝑓 + 𝐾 ∑︁ 𝑘 = { 𝑤 𝑣𝑐 (cid:107) v 𝑘 (cid:107) + (cid:107) 𝑊 𝑡𝑟 ( z 𝑘 − ¯ z 𝑘 )(cid:107) } (58)s.t. x 𝑘 + = 𝐴 𝑘 x 𝑘 + 𝐵 − 𝑘 u 𝑘 + 𝐵 + 𝑘 u 𝑘 + + Σ 𝑘 𝑇 + w 𝑘 + v 𝑘 ∀ 𝑘 = , . . . , 𝐾 − | x 𝑘 | < x 𝑚𝑎𝑥 ∀ 𝑘 = , . . . , 𝐾 (60) | u 𝑘 | < u 𝑚𝑎𝑥 ∀ 𝑘 = , . . . , 𝐾 (61) x = ¯ x 𝑖𝑛𝑖𝑡𝑖𝑎𝑙 (62) x 𝐾 = ¯ x 𝑓 𝑖𝑛𝑎𝑙 (63)A similar process is used to transcribe equations (17) - (24) into the Minimum-Effort Multi-Target Pointing OPT :minimize { x 𝑘 , u 𝑘 , v 𝑘 } 𝐾𝑘 = ∑︁ 𝑘 ∈ K (cid:8)(cid:13)(cid:13) ¯ q + 𝑘 q 𝑘 − q I (cid:13)(cid:13) + 𝛾 (cid:107) 𝝎 𝑘 − ¯ 𝝎 𝑘 (cid:107) (cid:9) + 𝐾 ∑︁ 𝑘 = (cid:40) 𝜌 (cid:107) u 𝑘 (cid:107) + 𝑤 𝑣𝑐 (cid:107) v 𝑘 (cid:107) + (cid:13)(cid:13)(cid:13)(cid:13)(cid:13) 𝑊 𝑡𝑟 (cid:32) (cid:34) x 𝑘 u 𝑘 (cid:35) − (cid:34) ¯ x 𝑘 ¯ u 𝑘 (cid:35) (cid:33)(cid:13)(cid:13)(cid:13)(cid:13)(cid:13) (cid:41) (64)s.t. x 𝑘 + = 𝐴 𝑘 x 𝑘 + 𝐵 − 𝑘 u 𝑘 + 𝐵 + 𝑘 u 𝑘 + + Σ 𝑘 𝑇 + w 𝑘 + v 𝑘 ∀ 𝑘 = , . . . , 𝐾 − | x 𝑘 | < x 𝑚𝑎𝑥 ∀ 𝑘 = , . . . , 𝐾 (66) | u 𝑘 | < u 𝑚𝑎𝑥 ∀ 𝑘 = , . . . , 𝐾 (67) x = ¯ x 𝑖𝑛𝑖𝑡𝑖𝑎𝑙 (68) x 𝐾 = ¯ x 𝑓 𝑖𝑛𝑎𝑙 (69)where we assume that the set of observation points is a subset of the discretization points, i.e., K ⊂ { , . . . , 𝐾 } .
5. Algorithm
The sequential convex programming algorithm used in this paper is listed in Algorithm 1 where S is an initial guessat the solution. As a notation convention used in the algorithm, superscript 𝑖 refers to the solution at the 𝑖 𝑡ℎ iteration ofthe algorithm: S 𝑖 : = { ¯ x 𝑖 , . . . , ¯ x 𝑖𝑘 , . . . , ¯ x 𝑖𝐾 } , { ¯ u 𝑖 , . . . , ¯ u 𝑖𝑘 , . . . , ¯ u 𝑖𝐾 } , ¯ 𝑡 𝑖𝑓 (70)The algorithm stops when either (1) the user-defined maximum number of SCP iterations 𝑁 𝑚𝑎𝑥 have been excuted, or(2) the algorithm has converged on a solution, where we define the convergence as the satisfaction of two conditions: ( 𝐽 𝑣𝑐 ≤ 𝜖 𝑣𝑐 ) ∧ ( 𝐽 𝑡𝑟 ≤ 𝜖 𝑡𝑟 ) (71)where 𝜖 𝑣𝑐 and 𝜖 𝑡𝑟 are user-defined convergence tolerances. The first condition ensures that a negligible amount ofvirtual controls is used, indicating that the converged solution is dynamically feasible. The second condition ensures9hat the solution remains sufficiently close to the nominal trajectory upon which the trajectory optimization problem wasformulated. Using a convergence flag, we keep track of whether the SCP subroutine has converged or not. Algorithm 1:
Sequential Convex Programming
Input : S Output : S 𝑖 , flag for 𝑖 = 𝑁 𝑚𝑎𝑥 do formulate OPT 𝑖 (cid:0) S 𝑖 − (cid:1) S 𝑖 ← OPT 𝑖 (cid:0) S 𝑖 − (cid:1) if S 𝑖 converged then flag = 1 returnif 𝑖 = 𝑁 𝑚𝑎𝑥 then flag = 0 return V. Numerical Examples
In this section we describe two examples to which we may apply the problem formulations stated in Section III. Wefirst demonstrate how the Minimum-Time Slew OCP is used to inform a constellation scheduler of the minimum timerequired by a reference satellite to conduct an arbitrary slew maneuver. We then assume that the constellation schedulerhas provided a desired target-pointing schedule and we plan an attitude reference trajectory using the Minimum-EffortMulti-Target Pointing OCP to realize the schedule. We model our reference satellite based on the physical parameterslisted in Table (1), representative of Planet’s Skysat, a satellite capable of agile maneuvering and imaging with sub-meter(50cm GSD) resolution [38].Parameter Value Units Description First Mention 𝑚
110 [kg] Satellite mass 𝑙 × 𝑤 × ℎ × ×
95 [cm] Satellite dimensions (cuboid) 𝐽 (cid:104) . . . . . . . . . (cid:105) [kg m ] Satellite inertia matrix Eqn (11) 𝑟 𝑚𝑎𝑥 𝑢 𝑚𝑎𝑥 𝐴 𝑟 (cid:104) -0 .
68 0 .
68 0 .
68 -0 . .
68 -0 .
68 0 .
68 0 . .
26 0 .
26 0 .
26 0 . (cid:105) [ ] Actuator Jacobian (4 rotors) Eqn (11) Table 1 Reference satellite parametersA. Minimum-Time Slew OCP Example
A target-sensing or data-downlinking schedule prescribes a sequence of desired instrument/antenna pointingdirections at specified times. To plan a schedule for each satellite in the constellation, the scheduler must considerthe time it takes for the satellite to slew from one orientation to another. For example, if the distance between twoconsecutively scheduled targets on the Earth’s surface is large, it may not be feasible for the satellite to slew from onepointing orientation to another given the time interval between points. The slew time depends not only on a satellite’sorbital motion and the relative distance between targets, but also on its mass distribution and actuator constraints. Toour knowledge, a closed-form expression for computing minimum slew times between any arbitrary orientations doesnot exist (even for a symmetric rigid body with independent three-axis control).Our approach is to apply the Minimum-Time Slew OCP (Section III.A) over a gridded space of 3D rotations. Since anyattitude trajectory between any two orientations can be summarized with a single rotation about some axis, we parameterizethe rotations using the Euler axis-angle representation. As illustrated in Fig. (1), we consider 100 equidistributed axesof rotation and 88 rotations of increasing magnitude: 𝜃 ∈ { -180 , -175 , . . . , -10 , -9 , -1 , , , . . . , , , . . . , } degreesto cover a relatively fine grid of the entire space of 3D rotations, resulting in 8,800 unique Euler axis-angle tuples.10 ig. 1 (Left) 100 equidistributed points on unit sphere, (Right) 8,800 rotations of a pointing vector. We set up 8,800 instances of Minimum-Time Slew OCP, where we set the initial value of the quaternion in (15) tobe the identity quaternion (i.e., q = q I ) and the final desired quaternion in (16) to be the quaternion representation ofan Euler axis-angle rotation. Using the sequential convex programming method, we efficiently solve each probleminstance off-line with a computation time on the order of seconds. More efficient SCP implementations [33] may solveeach problem in milliseconds and the overall procedure can be parallelized for even faster off-line computation. Byrecording the minimum maneuver time and ADCS energy consumption from the solution to each problem instance, wecan produce lookup tables or data-fitted functions to be used by a scheduler. B. Minimum-Effort Multi-Target Pointing OCP Example
Given a desired attitude pointing schedule { 𝑡 𝑘 , ¯ q ( 𝑡 𝑘 ) , ¯ 𝝎 ( 𝑡 𝑘 )} ∀ 𝑘 ∈ K , we may plan a continuous trajectory thatpasses through each attitude point in the sequence by using Minimum-Effort Multi-Target Pointing OCP (Section III.B).To motivate this method, we consider the remote-sensing application described in [2], where a 24-satellite, 3-planeWalker-Delta constellation is simulated at 710 km altitude, 98.5 deg inclined, circular orbits over a 6-hour durationwith the orbital mechanics module of the D-SHIELD software suite [4]. Simulation results include the orbital statesof the satellites at each time step as well as access times when user-defined target grid points (described by lattitudeand longitude coordinates) are observable by a satellite. In our example we consider 42 urban regions that experiencefrequent episodic precipitation and are prone to flooding [2]. Each of these globally-distributed watersheds covers an 80km area spanned by 121 grid points. Fig. 2 (Left) Satellite groundtrack and 16 observable watersheds, (Right) Cluster of 6 watersheds
In our example, we focus on a single satellite of the constellation and observe in the first plot of Fig. (2) that thesatellite can access 16 of the 42 watersheds (each containing 121 green grid points) during the 6-hour time span. Onlya subset of the 42 watersheds are accessible due to access restrictions, bounds on maximum off-pointing angles, orocclusion by the Earth’s surface. The second plot of Fig. (2) provides a zoomed-in view of the groundtrack, containinga dense cluster of 6 regions with a total of 726 potentially accessible target grid points. As the satellite passes over this11ense cluster, the satellite may be commanded by the scheduler to perform a rapid sequence of agile slewing maneuversto acquire, point, and track desired target grid points. The commanded schedule during this part of the orbital trajectorymay be the most difficult to execute by an attitude control system, requiring high-frequency intra-region slewing as wellas large-angle, inter-region slewing. Hence, our example focuses on fulfilling a pointing schedule that maximizes thenumber of observed grid points in this 6-region cluster, treating it as a stress test for our proposed approach to executingmulti-target pointing schedules.Figure (3) shows both Earth-Centered Inertial (ECI) and Earth-Centered Earth-Fixed (ECEF) views of the satellite’sorbit, where the red vectors point from the satellite to green grid points on the Earth’s surface. We note that our stresstest focuses on a subset (i.e., a dense 6-region cluster) of the 16 accessible regions shown in Fig. (3) . At any giventime, a desired pointing frame may be uniquely defined with its z-axis aligned with a red pointing vector, x-axis inthe orbital plane towards the direction of motion, and y-axis completing the right-hand triad. The orientation of thisframe with respect to the ECI frame is the desired pointing attitude to acquire. When a desired attitude is defined ata particular time instance, we may also define the desired body angular velocity to be zero for an interval about thespecified time instance (e.g., if ∃ ¯ q ( 𝑡 𝑘 ) , then { ¯ 𝝎 ( 𝑡 𝑘 − ) , ¯ 𝝎 ( 𝑡 𝑘 ) , ¯ 𝝎 ( 𝑡 𝑘 + )} = ). Specifying zero velocity at the time oftarget observation may mitigate the effects of motion blur while imaging. Fig. 3 (Left) Earth-Centered Inertial, (Right) Earth-Centered Earth-Fixed coordinate systems
Using the following guidelines, we manually design a sequence of desired pointing orientations (i.e., an attitudepointing schedule) to maximize grid point (GP) observations across the dense 6-region cluster: • Include as many of the 726 GPs as feasibly possible given minimum attitude slew time between observations. • Visit each of the six regions (or as many as possible). • Finish observing a region before moving to the next region (i.e., cannot return to a previously observed regions). • Observe GPs using a “sweeping” pattern, similar to the whiskbroom pattern used by early Landsat satellites [39]. • The first accessible points in a region are observed first, deciding direction of the intra-region, sweeping pattern. • The first accessible regions are observed first, deciding direction of the inter-region slewing.As the satellite passes over the 6-region cluster, it has an approximately 10-minute window from the time when oneof the 726 GPs first becomes accessible to the time that none are accessible. From an altitude of 710km, we estimatethat the satellite requires at least two seconds to conduct a minimum-time, rest-to-rest slew from a grid point in thenadir direction to an adjacent grid point 8 km away. We also estimate that it can require at least 15 seconds to slew theseveral hundreds of kilometers between certain pairs of regions in the cluster. These estimates for slew time betweenneighboring GPs (intra-region slewing) and neighboring regions (inter-region slewing) are based on insights gainedfrom Section (V.A). We use the slew-time estimates to design a pointing schedule with a one-second sample time.The first plot of Fig. (4) shows a manually-designed schedule that visits five of the six regions. In the second plotof Fig. (4), we illustrate the “sweeping” pattern that we adopt within each region to observe 33 GPs per region. Ared marker signifies a GP chosen for observation and a green marker is an accessible GP not chosen for observation.The red path illustrates the sequence in which the GPs are observed as well as the approximate trajectory of thesatellite’s pointing vector on the Earth’s surface. With 33 GPs in five regions, we design a schedule that covers atotal of 165 GPs or 22.7 percent of the total number of accessible GPs. In our schedule, we omit the region locatednear 25 ◦ lattitude and 85 ◦ longitude since its inclusion reduces the total number of GP observations, assuming the12 ig. 4 (Left) Inter-region slewing, (Right) Intra-region slewing between 33 red GPs selected per region. same number of GP observations per region. Given the 10-minute window and the minimum-slew time requirementsbetween GPs and between regions, we find that the 165-GP pointing schedule is ambitious and close to the maximumnumber of grid points that a satellite may feasibly observe across the 6-region cluster. The other grid points mustbe observed on a subsequent revisit or by another satellite. Compared to this heuristics-based schedule, we expectthat an optimization-based schedule could further maximize the number of GP observations. An optimization-basedscheduler may explicitly use the minimum maneuver time model described in Section V.A to accurately estimate thefastest transition time between GP observations. By finding the most time-efficient sequence of GPs, the scheduler canmaximize the number of observations or maximize science utility [4]. VI. Results
We apply the trajectory optimization method described Section IV to the examples described in Sections V.A andV.B. For each problem, we start the SCP Algorithm (1) with a crude initial guess S , as represented in (70), that consistsof a spherical linear interpolated [40] trajectory between desired quaternions, zero body angular velocity, zero rotormomenta, and zero rotor torque. We recall that our reference satellite is modeled after Planet’s Skysat with representativephysical parameters listed in Table (1). In Table (2), we list the parameters used in each problem formulation and notethat the Minimum-Time Slew OPT and Minimum-Effort Multi-Target Pointing OPT use 𝐾 =
30 and 601 discretizationnodes (38), respectively. We choose a relatively small number of nodes for the Minimum-Time Slew OPT since weexpect that it is sufficient for the short time scales in this class of problems. That is, a minimum-time slew between anytwo orientations should be on the order of seconds. For the Minimum-Effort Multi-Target Pointing OPT, we specify afixed final time of 600 seconds and use 601 temporal nodes so that each node represents one second of the ten-minuteschedule. We then add cost terms (64) for each node that corresponds to a time when a desired attitude is specified.Problem Formulation SCP Algorithm Convergence ResultsOPT
𝐾 𝛾 𝜌 𝑁 𝑚𝑎𝑥 𝑤 𝑣𝑐 𝑤 𝑡𝑟 𝜖 𝑣𝑐 𝜖 𝑡𝑟 𝑁 𝑐𝑜𝑛𝑣𝑒𝑟𝑔𝑒 𝑡 𝑐𝑜𝑛𝑣𝑒𝑟𝑔𝑒 Min-Time ( 𝜃 = ◦ ) 30 - - 20 1e5 1e-1 1e-5 1e-5 10 20.3 secMin-Error, Mult-Pt 601 1e5 1e0 20 1e5 1e-1 1e-3 1e-4 9 144.6 sec(Eqn No.) (38) (64) (64) Algo (1) (51) (53) (71) (71) - - Table 2 SCP Algorithm and Objective Function parameters, Convergence results
For both OPTs, the SCP algorithm reaches convergence, as defined by equations (51) and (52), the convergencecondition (71), and the SCP parameters listed in Table 2. Due in part to the small number of nodes used in theMinimum-Time Slew OPT, each of the 8,800 problem instances are solved in a matter of seconds. However, theMinimum-Effort Multi-Target Pointing OPT, with 601 nodes, requires 144.6 seconds to converge to a solution. Figure(5) illustrates how the virtual control and trust region penalty terms for the Minimum-Effort Multi-Target PointingOPT decrease below user-defined tolerance values by the 9th SCP iteration. The small value for 𝐽 𝑣𝑐 signifies that thedynamics, state constraints and input constraints are satisfied using the input trajectory { u ( 𝑘 )} 𝐾𝑘 = and only a negligible13ontribution from virtual controls (Section IV.D.1). Furthermore, the small value for 𝐽 𝑡𝑟 shows that the solutionat the 9th iteration has not changed significantly from that of the previous iteration. This signifies that the currentsolution is based on an accurate linearization of the nonlinear dynamics about the previous solution, and that negligibleimprovement may be found in further iterations. The problems are written in MATLAB and solved using CVX, apackage for specifying and solving convex programs [41], on a 2015 MacBook Pro with 2.5 GHz Intel Core i7 processor. Iteration Number
Jvc -3 J t r Fig. 5 Virtual Control and Trust Region ( 𝐽 𝑣𝑐 , 𝐽 𝑡𝑟 ) penalty terms fall under tolerance values ( 𝜖 𝑣𝑐 , 𝜖 𝑡𝑟 ) by the9th iteration, representing SCP algorithm convergence to a solution.A. Minimum-Time Slew OPT Result We first study the results of applying the Minimum-Time Slew OPT to a particular slew maneuver: a 60-degreerotation about an arbitrary non-principal axis. In Fig. (6), we illustrate the spacecraft’s motion where the blue vectorrepresents the body z-axis, aligned with an instrument’s pointing vector. The trajectories traced out by the body axesrepresent the minimum-time slew maneuver to achieve the reorientation. We qualitatively note that this maneuver aboutan arbitrary non-principal axis does not resemble an eigenaxis slew. However, for slew maneuvers about any of theprincipal axes, we confirm that the minimum-time maneuver is an eigenaxis slew aligned with the principal axis.
Fig. 6 Minimum-time slew maneuver (60-degree rotation about an arbitrary axis)
We characterize the solution for this particular problem instance using five quantities: maneuver time, rotor momenta,rotor torque, ADCS instantaneous power drawn and ADCS cumulative energy consumption. The time series of thesequantities are shown below. In the first plot of Fig. (7), we note that the quaternion trajectory starts from a nominalorientation (represented by the identity quaternion) and reaches the desired terminal nodes, corresponding to thegiven axis-and-angle rotation. The second plot of Fig. (7) shows that the body starts and ends with zero angularvelocity, representing a rest-to-rest slew maneuver. The angular speed is greatest about the body z-axis, confirmingthe predominantly yaw-axis rotation depicted in Fig. (6). Although we have not enforced bounds on maximum body14ngular speed (i.e., body slew rate) in our results, we may easily add these constraints to the problem formulation. Forexample, to maintain star tracker accuracy we may impose an upper bound on the body slew rate.The optimal maneuver time is approximately 21.3 seconds. Given the same axis-and-angle and actuator constraints,we expect that a spacecraft with less mass but proportionate mass distribution requires less maneuver time. Hence, thetrend towards small satellites not only results in cost efficiency but also more nimble operation. t [s] q [ ] q1q2q3q4 t [s] -20246 w [ deg / s ] wxwywz Fig. 7 (Left) Quaternion parameterization of attitude, (Right) Body Angular Velocity
The first plot of Figure (8) shows that the rotor momenta do not saturate for this 60-degree maneuver. For largerrotation magnitudes, we expect that the momentum trajectories may approach the magnitude bound of 0.80 Nms,impacting maneuver time. We note the non-zero rotor momenta at the end of the maneuver caused by the input torques.In practice, momentum will also build up in the rotors as they counteract torques due to environmental disturbances.This momentum buildup may be dumped with magnetorquers or other compensating mechanism. In the second plot,we observe that the rotor torque magnitudes saturate at 0.06 Nm and follow a “bang-bang” structure, consistent withclassical minimum-time maneuvers [42]-[43]. This observation confirms our intuition that performing a rotation inminimum time requires the actuators to perform at their limits. t [s] -1-0.500.51 r [ N m s ] rotor 1rotor 2rotor 3rotor 4rmax t [s] -0.06-0.04-0.0200.020.04 u [ N m ] rotor 1rotor 2rotor 3rotor 4umax Fig. 8 (Left) Rotor Momentum, (Right) Rotor Torque
Finally, in Fig. (9), we note that the peak instantaneous power drawn by the ADCS stays below 15 W and that thetotal energy consumption is less than 150 J for the maneuver. Given mission constraints on the instantaneous power orenergy, we may enforce constraints (27) and (28). t [s] P [ W ] rotor 1rotor 2rotor 3rotor 4total t [s] E [ J ] rotor 1rotor 2rotor 3rotor 4total Fig. 9 (Left) Instantaneous Power drawn, (Right) Cumulative Energy Consumed
The results above are for one of the 8,800 problem instances described in (V.A). After solving each probleminstance, we plot the corresponding slew maneuver time and ADCS energy consumption as functions of the rotationmagnitude in Fig. (10). For each of the 100 equidistributed axes of rotation, we observe that optimal maneuver time is apower function of the rotation magnitude and that required energy is a linear function of rotation magnitude. Theserelationships allow us to approximate the minimum maneuver time and required energy consumption to conduct anyarbitrary rotation by the reference satellite described in Table (1). The relationships may be represented in lookup tablesor function approximators that can be quickly evaluated by an on-line constellation scheduler.15
200 -100 0 100 200
Rotation Magnitude [deg] M anue v e r T i m e [ s ] Body x-axis fitBody y-axis fitBody z-axis fit -200 -100 0 100 200
Rotation Magnitude [deg] R equ i r ed E ne r g y E [ J ] Body x-axis fitBody y-axis fitBody z-axis fit
Fig. 10 (Left) Minimum-time, (Right) Energy as functions of eigenaxis rotation magnitude 𝜃 To empirically verify the optimality of the solutions found by our approach, we solve the Minimum-Time Slew OPTfor rotations about the principal axes (i.e., about body x-, y-, and z-axes). Using least-squares fitting, we determinethe coefficients for a minimum-time power model, 𝑇 = 𝑎𝜃 𝑏 , and an energy linear model, 𝐸 = 𝑐𝜃 + 𝑑 , valid foreigenaxis rotations | 𝜃 | > ◦ about each axis. As shown in Table (3), the power coefficient in the minimum-time modelapproximately takes on the value of 𝑏 = . Table 3 Coefficients produced by least squares fitting (valid for 𝜃 > ◦ ) Beyond double-integrator dynamics about principal axes, our approach can be applied to high-fidelity spacecraftmodels and arbitrary rotations. Since the Minimum-Time Slew OPT can be applied to non-rest-to-rest and non-principalaxis rotations, where nonlinear dynamical coupling exists, the minimum-slew-time/energy estimates found by ourapproach may be more accurate than those based on the classical minimum-time, rest-to-rest eigenaxis slew solution.Furthermore, we may explicitly consider system- or actuator-specific constraints, such as slew rate or rotor momentumbounds, that are not included in the classical formulation.
B. Minimum-Effort Multi-Target Pointing OPT Result
We apply the Minimum-Effort Multi-Target Pointing OPT to the 165-GP attitude pointing schedule for the 6-regioncluster example described in Section (V.A). We recall that the schedule is a discrete sequence of desired pointingorientations: { 𝑡 𝑘 , ¯ q ( 𝑡 𝑘 ) , ¯ 𝝎 ( 𝑡 𝑘 )} ∀ 𝑘 ∈ K , where | K | = t [min] -1-0.500.51 q [ ] q1q2q3q4 Fig. 11 Quaternion parameterization of attitude trajectory and the desired quaternion schedule it attempts to satisfy. Equation (73) determines the average error acrossall points. As recorded in Table 5, we compute the maximum point-wise quaternion error to be 𝑞 𝑚𝑎𝑥𝑒 = . 𝑞 𝑎𝑣𝑔𝑒 = . 𝑞 𝑚𝑎𝑥𝑒 = max 𝑘 ∈ K (cid:13)(cid:13) ¯ q + 𝑘 q 𝑘 − q I (cid:13)(cid:13) (72) 𝑞 𝑎𝑣𝑔𝑒 = | K | ∑︁ 𝑘 ∈ K (cid:13)(cid:13) ¯ q + 𝑘 q 𝑘 − q I (cid:13)(cid:13) (73)In Fig. (12), we show the corresponding angular velocity trajectory where we observe that the spacecraft reachesnearly-zero angular velocity at all 165 desired points. Similar to the performance metrics for the quaternion trajectory,we also introduce equations (74) and (75) to quantify the maximum point-wise error and the average error in the angularvelocity trajectory with respect to the desired schedule. We record a maximum point-wise error of 𝜔 𝑚𝑎𝑥𝑒 = . 𝜔 𝑎𝑣𝑔𝑒 = . 𝜔 𝑚𝑎𝑥𝑒 = max 𝑘 ∈ K (cid:107) 𝝎 𝑘 − ¯ 𝝎 𝑘 (cid:107) (74) 𝜔 𝑎𝑣𝑔𝑒 = | K | ∑︁ 𝑘 ∈ K (cid:107) 𝝎 𝑘 − ¯ 𝝎 𝑘 (cid:107) (75)To further reduce the angular velocity error, we may increase the relative state penalty 𝛾 in the Minimum-EffortMulti-Target Pointing OPT objective (64). We may also choose to include an explicit constraint on the error using (26)and attempt to find a feasible solution.We also note in Fig. (12) that the angular velocity components peak while slewing between the five regions. Thelarge peaks are due to the limited time spans allotted for the inter-region slews specified in our manually-designedschedule. Imposing body slew rate constraints would reduce the amplitude of the peaks but may also result in increasedquaternion error, especially for the initial points of each region when high slew rates are required to move quickly fromthe previous region. A “relaxed” schedule that allows for more time to conduct inter-region slews could satisfy slew rateconstraints with minimal quaternion error but at the expense of observing more GPs.Quaternion Error Angular Velocity Error 𝑞 𝑚𝑎𝑥𝑒 𝑞 𝑎𝑣𝑔𝑒 𝜔 𝑚𝑎𝑥𝑒 𝜔 𝑎𝑣𝑔𝑒 Table 4 Reference trajectory error with respect to desired GP schedule
As shown in Fig. (13), the rotor momenta stays within the 0.8 Nms bounds and peaks when conducting large anglemaneuvers, such as slewing between regions. Despite the large distances between certain pairs of regions, the rotormomenta do not saturate. This suggests that the reference spacecraft is capable of slewing between regions with muchlarger ground separations before the rotor momenta saturate.17 t [min] -6-4-20246 w [ deg / s ] wxwywz Fig. 12 Body Angular Velocity t [min] -0.8-0.6-0.4-0.200.20.40.60.8 r [ N m s ] rotor 1rotor 2rotor 3rotor 4rmax Fig. 13 Rotor Momentum
Despite the control effort penalty in the objective (64), which encourages minimal control effort, we observe that therotor torque trajectories saturate at the 0.06 Nm bounds in Fig. (14) in a “bang-bang,” minimum-time fashion. Thissuggests that the spacecraft is working at its performance limits to achieve an aggressive pointing schedule. We expectactuator saturation to occur for schedules that maximize the number of observed grid points in a certain amount of time. u [ N m ] u u t [min] -0.0500.05 u Fig. 14 Rotor Torque
The rotor momentum and torque trajectories inform us of the reference spacecraft’s agile performance limits.Schedules requiring large-angle slew maneuvers may saturate rotor momentum bounds, limiting how far the spacecraft’spointing vector can span in a given amount of time. Schedules requiring rapid slewing between attitude points maysaturate rotor torque bounds. Roughly speaking, for a given spacecraft design (i.e., mass moment-of-inertia, actuatorconfiguration), rotor momentum and torque bounds dictate spatial range-of-motion and speed, respectively.Similar to the body angular velocity and rotor momentum trajectories, we observe in Fig. (15) that the totalinstantaneous power drawn by the rotors (i.e., ADCS) peaks when slewing between regions. We recall that if maximumpower bounds must be enforced on the ADCS, we may add constraint (27) in our problem formulation.18 t [min] P [ W ] rotor 1rotor 2rotor 3rotor 4total Fig. 15 ADCS Instantaneous Power drawn
In Fig. 16, we note that the ADCS requires approximately 550 J to conduct this ten-minute, minimum-effortmulti-target pointing trajectory. If this trajectory is too energetically expensive, we may increase the relative controleffort weight 𝜌 in the problem objective (17) to further penalize energy consumption. However, the choice may come atthe expense of increased attitude error. If a specific energy bound must be satisfied, we may explicitly enforce it byadding constraint (28). In the case that a feasible trajectory with an energy upper bound cannot be found, the attitudepointing schedule may need to be relaxed by observing few regions and/or GPs. t [min] E [ J ] rotor 1rotor 2rotor 3rotor 4total Fig. 16 ADCS Cumulative Energy consumedC. Open-loop versus Closed-loop Performance
Using the solution to the
Minimum-Effort Multi-Target Pointing OPT , we can apply the reference torque trajectory { u 𝑂𝐿 𝑘 : = ¯ u 𝑘 } 𝐾𝑘 = as a linearly-interpolated, open-loop control command in simulation with the dynamics described in (3)and parameters listed in Table (1). We observe that the resulting state trajectory matches the solution { ¯ q 𝑘 , ¯ 𝝎 𝑘 , ¯ r 𝑘 } 𝐾𝑘 = with negligible error, signifying that the OPT approximation, solved using sequential convex programming, has produceda feasible solution for the original
OCP formulation. However, in the presence of external disturbances, parameteruncertainty or unmodeled dynamics, a feedback control law must be applied to track the desired attitude trajectory.We design an optimal tracking controller based on the solution to the LQR problem (see Appendix) and use theresulting control law as feedback term to regulate deviations from the desired trajectory. The rotor torque command thatwe apply to (3) consists of both the feedforward and feedback terms: { u 𝐶𝐿 𝑘 : = ¯ u 𝑘 + 𝛿 u 𝑘 } 𝐾 − 𝑘 = , where 𝛿 u is a time-varyingstate feedback law.To compare open-loop versus closed-loop performance in presence of model mismatch, we perturb the inertia matrixused in simulation. While we assume a nominal inertia matrix 𝐽 in the OCP formulation and the LQR tracking controllaw design, we use a significantly perturbed ˜ 𝐽 in simulation: 𝐽 = . . . . . . . . . kg · m ˜ 𝐽 = . . . . . . . . . kg · m (76)19s illustrated in Fig. 17, applying the optimal torque reference trajectory in open-loop produces the dashedquaternion trajectory, which fails to pass through the desired quaternion points specified by our schedule. However,when using the LQR tracking controller, the closed-loop trajectory is able to pass through most of the desired quaternionpoints. In Fig. 18 we confirm that the amplitude of the angular velocity peaks during inter-region slewing are reduced. t [min] -1-0.500.51 q [ ] desiredopen-loopLQR tracker Fig. 17 LQR closed-loop and open-loop quaternion trajectories t [min] -8-6-4-202468 w [ deg / s ] desiredopen-loopLQR tracker Fig. 18 LQR closed-loop and open-loop angular velocity trajectories
In Section VI.B, we use the performance metrics (72)-(75) to measure the error between the optimized referencetrajectory and the desired attitude schedule. In this section, we use the same metrics to compute the error of theclosed-loop (or open-loop) simulated trajectory with respect to the reference trajectory. As shown in Table 5, theQuaternion Error Angular Velocity ErrorControl Strategy 𝑞 𝑚𝑎𝑥𝑒 𝑞 𝑎𝑣𝑔𝑒 𝜔 𝑚𝑎𝑥𝑒 𝜔 𝑎𝑣𝑔𝑒 Open-loop 0.5947 0.3417 3.5724 0.4116LQR tracker 0.0653 0.0085 2.0227 0.3515
Table 5 Simulation trajectory error with respect to reference trajectory
LQR-based, closed-loop trajectories result in smaller error across all four performance metrics. In particular, themaximum point-wise quaternion (i.e., attitude) error is reduced by an order of magnitude when using the trackingcontroller. In Fig. 19, we plot the quaternion and angular velocity errors as functions of the simulation time. We observethat the maximum point-wise error in the open-loop quaternion trajectory coincides with the final point in the trajectory.Since the error accumulates over time, we can infer that schedules prescribed for long durations of time will result inlarge attitude errors if optimal torque commands are applied in open-loop. In contrast, the LQR tracker significantlyreduces error in the closed-loop quaternion trajectory. Furthermore, both the maximum and average angular velocityerror is reduced with the LQR tracker. Further reduction in the angular velocity error may be achieved (at the expense ofeither increased quaternion error or control effort) by tuning the LQR weight matrices (80).20 q e [ ] open-loopLQR tracker t [min] w e [ deg / s ] Fig. 19 (Top) Quaternion error, (Bottom) Angular velocity error
VII. Conclusion
We have formulated two optimal control problems, the Minimum-Time Slew OCP and the Minimum-EffortMulti-Target Pointing OCP, that can be applied to high-fidelity spacecraft models with practical mission objectivesand constraints. The 3-DOF, rotational motion of a spacecraft may be modeled with nonlinear gyrostat and actuatordynamics, quaternion kinematics, arbitrary physical parameters, and arbitrary actuator configurations. Constraintssuch as bounds on rotor torque, rotor momentum, body slew rate, ADCS power and energy are included. Usingstate-of-the-art techniques in sequential convex programming, we transcribe the non-convex, continuous-time optimalcontrol problems into convex, finite-dimensional optimization problems (OCP → OPT) to be solved efficiently.The Minimum-Time Slew OPT accurately estimates the required time and energy to conduct time-optimal, rest-to-restmaneuvers between any two arbitrary orientations. Compared with simulation-based approaches that implementnear-optimal, eigenaxis-slew control laws about body-fixed control axes, our trajectory optimization approach is notrestricted to eigenaxis slews and can be applied to more general spacecraft models and constraints. The time and energyestimates can be computed on a grid of the space of 3D rotations to produce models for use by a constellation scheduler.Given a desired schedule of discrete pointing orientations, the Minimum-Effort Multi-Target Pointing OPT producesa dynamically-feasible, continuous reference trajectory that can achieve the pointing schedule with minimal effort. Tomitigate the effects of model mismatch or external disturbances, a trajectory-tracking controller is used to regulatedeviations from the trajectory. Compared to existing approaches that step from one discrete reference set-point toanother, our method allows us to design the entire multi-target pointing trajectory as a single unabridged maneuver.We have applied both problems formulations to an example multi-target observation mission and reference spacecraft,demonstrating their potential use by a remote sensing constellation scheduler.The two problem formulations are complementary and may significantly elevate the performance of agile, remotesensing satellite constellations. With more accurate estimates of the minimum time and energy required to conductparticular reorientations, a constellation scheduler is better informed to design aggressive multi-target pointing schedulesthat maximize constellation utility (e.g., quantity or science value of observations). The satellites of the constellationsmay then execute the schedules by planning feasible attitude trajectories that explicitly consider state and input constraints.Furthermore, given recent advances in the computing power of small satellite on-board computers, real-time solverimplementations allow a constellation to respond quickly to changes in the environment or mission objectives.Future work will include a rigorous comparison of our proposed method to existing approaches in the presence ofmodel uncertainty, disturbances, and the use of state estimation. We will also study other applications beyond groundtarget observation, including the tracking of cloudbows and aerosols in the atmosphere. Finally, our attitude trajectoryplanning approach will be integrated with an optimization-based scheduler so we may compare its performance to themanually-designed, baseline schedule used in this paper.
Appendix : Linear Quadratic Trajectory Tracking Controller Design
We use the closed-form solution of the (unconstrained) finite-horizon, discrete-time LQR problem [42] - [44] todesign a time-varying, state feedback law that regulates deviation from desired trajectory in closed-loop simulation.Due to the intrinsic unit-norm constraint on the quaternion, we convert the error-quaternion to the Euler axis and21ngle representation with a transformation 𝝓 = ℎ ( ¯ q + q ) in our tracking controller implementation. Details on thetransformation may be found in [34]-[36].Using error variables 𝝓 , 𝛿 𝝎 : = ( 𝝎 − ¯ 𝝎 ) and 𝛿 r : = ( r − ¯ r ) to represent deviations from the desired trajectory: 𝛿 x ( 𝑡 ) : = [ 𝝓 ( 𝑡 ) (cid:62) 𝛿 𝝎 ( 𝑡 ) (cid:62) 𝛿 r ( 𝑡 ) (cid:62) ] (cid:62) , 𝛿 u ( 𝑡 ) : = u ( 𝑡 ) − ¯ u ( 𝑡 ) , (77)the error system dynamics are: (cid:164) 𝛿 x ( 𝑡 ) = ˆ 𝐴 ( 𝑡 ) 𝛿 x ( 𝑡 ) + ˆ 𝐵 ( 𝑡 ) 𝛿 u ( 𝑡 ) , (78)ˆ 𝐴 ( 𝑡 ) : = −[ ¯ 𝝎 ( 𝑡 )] × I − 𝐽 − (cid:0) [ ¯ 𝝎 ( 𝑡 )] × 𝐽 − [ 𝐽 ¯ 𝝎 ( 𝑡 ) + 𝐴 𝑟 ¯ r ( 𝑡 )] × (cid:1) − 𝐽 − [ ¯ 𝝎 ( 𝑡 )] × 𝐴 𝑟 , ˆ 𝐵 ( 𝑡 ) : = − 𝐽 − 𝐴 𝑟 I , (79)where and I are the 3 × [ ] × described in (4). In the model above, note that we have approximated the error-quaternion kinematics using the Euleraxis and angle representation. Assuming a sample time of Δ 𝑡 = 𝑡 𝑓 𝐾 − , we discretize the error dynamics with a zero-orderhold on 𝛿 u . The attitude trajectory tracking problem is formulated as:minimize { 𝛿 x 𝑘 , 𝛿 u 𝑘 } 𝐾𝑘 = 𝛿 x (cid:62) 𝐾 𝑄 𝐾 𝛿 x 𝐾 + 𝐾 − ∑︁ 𝑘 = (cid:8) 𝛿 x (cid:62) 𝑘 𝑄 𝑘 𝛿 x 𝑘 + 𝛿 u (cid:62) 𝑘 𝑅𝛿 u 𝑘 (cid:9) (80)s.t. 𝛿 x 𝑘 + = ˆ 𝐴 𝑘 𝛿 x 𝑘 + ˆ 𝐵 𝑘 𝛿 u 𝑘 ∀ 𝑘 = , . . . , 𝐾 − 𝛼 : 𝑄 𝑘 : = (cid:40) 𝑄, 𝑘 ∉ K 𝛼 · 𝑄, 𝑘 ∈ K where 𝑄 (cid:60) , 𝑅 (cid:31) , 𝛼 > . (82)The problem formulation admits a unique closed-form solution with feedback control law { 𝛿 u 𝑘 : = − K 𝑘 𝛿 x 𝑘 } 𝐾 − 𝑘 = , where: K 𝑘 = ( 𝑅 𝑘 + ˆ 𝐵 (cid:62) 𝑘 𝑃 𝑘 + 𝐵 𝑘 ) -1 ˆ 𝐵 (cid:62) 𝑘 𝑃 𝑘 + ˆ 𝐴 𝑘 (83)and 𝑃 𝑘 is found recursively from 𝑃 𝐾 = 𝑄 𝐾 using the following discrete-time dynamic Riccati equation: 𝑃 𝑘 = K (cid:62) 𝑘 𝑅 𝑘 K 𝑘 + ( ˆ 𝐴 𝑘 − ˆ 𝐵 𝑘 K 𝑘 ) (cid:62) 𝑃 𝑘 + ( ˆ 𝐴 𝑘 − ˆ 𝐵 𝑘 K 𝑘 ) + 𝑄 𝑘 . (84) Acknowledgments
This project has been funded and supported by NASA’s Earth Science Technology Office via the AdvancedInformation Systems Technology grant. We thank the Orbital R&D Team at Planet Labs for providing representativeparameter values for the reference satellite used in this paper and for sharing insight on the state-of-the-practice insatellite attitude control.
References [1] S. Nag, A. Li, J.H. Merrick, “Scheduling algorithms for rapid imaging using agile Cubesat constellations,” in
Advances inSpace Research , vol. 61, no. 3, pp. 891-913, November 2017.[2] S. Nag, A. Li, V. Ravindra, M.S. Net, K.M. Cheung, R. Lammers, B. Bledsoe, “Autonomous Scheduling of Agile SpacecraftConstellations with Delay Tolerant Networking for Reactive Imaging,” in
Proceedings of the International Conference onPlanning and Scheduling SPARK Workshop , Berkeley, CA, USA, July, 2019.[3] S. Nag, M.S. Net, A. Li, V. Ravindra, “Designing a Disruption Tolerant Network for Reactive Spacecraft Constellations,” in
Proceedings of the ASCEND , Virtual Event, November 2020.[4] S. Nag, A. Aguilar, R. Akbar, A. Azemati, J. Frank, R. Levinson, A. Li, M. Moghaddam, V. Ravindra, D. Selva, “D-SHIELD:Distributed Spacecraft with Heuristic Intelligence to Enable Logistical Decisions,” in
Proceedings of the IEEE InternationalGeoscience and Remote Sensing Symposium , Virtual Event, September 2020.
5] V. Shah, V. Vittaldev, L. Stepan, C. Foster, “Scheduling the World’s Largest Earth-Observing Fleet of Medium-ResolutionImaging Satellites,” in
Proceedings of the International Workshops on Planning and Scheduling for Space , Berkeley, CA, USA,July, 2019.[6] M. Harris, "Tech giants race to build orbital internet [News]," in
IEEE Spectrum , vol. 55, no. 6, pp. 10-11, June 2018.[7] A. Li, J. Mason, “Optimal Utility of Satellite Constellation Separation with Differential Drag,” in
Proceedings of the AIAASpace Forum , San Diego, CA, USA, August 2014.[8] C. Foster, J. Mason, V. Vittaldev, L. Leung, V. Beukelaers, L. Stephan, R. Zimmerman, “Constellation Phasing with DifferentialDrag on Planet Labs Satellites,”
Journal of Spacecraft and Rockets , vol. 55, no. 2, pp. 473-483, March 2018.[9] E. Sin, M. Arcak, A. Packard, “Small Satellite Constellation Separation using Linear Programming based Differential DragCommands,” in
Proceedings of the IEEE American Control Conference , Milwaukee, WI, USA, June 2018.[10] E. Sin, H. Yin, M. Arcak, “Passivity-based Distributed Acquisition and Station-keeping Control of a Satellite Constellation inAreostationary Orbit,” in
Proceedings of the ASME Dynamic Systems and Control Conference , Virtual Event, October 2020.[11] B. Wie, H. Weiss, A. Arapostathis, “Quaternion Feedback Regulator for Spacecraft Eigenaxis Rotations,”
Journal of Guidance,Control, and Dynamics , vol. 12, no. 3, pp. 375-380, May 1989.[12] J.T. Wen, K. Kreutz-Delgado, “The Attitude Control Problem,”
IEEE Transactions on Automatic Control , vol. 36, no. 10, pp.1148-1162, October 1991.[13] K.D. Bilimoria, B. Wie, “Time-Optimal Three-Axis Reorientation of a Rigid Spacecraft,”
Journal of Guidance, Control, andDynamics , vol. 16, no. 3, pp. 446-452, May 1993.[14] B. Wie, J. Lu, “Feedback Control Logic for Spacecraft Eigenaxis Rotations Under Slew Rate and Control Constraints,”
Journalof Guidance, Control, and Dynamics , vol. 18, no. 6, pp. 1372-1379, November 1995.[15] B. Wie, D. Bailey, C. Heiberg, “Rapid Multi-Target Acquisition and Pointing Control of Agile Spacecraft,” in
Proceedings ofthe AIAA Guidance, Navigation, and Control Conference , Denver, CO, USA, August 2000.[16] J.T. Betts, “Survey of Numerical Methods for Trajectory Optimization,”
Journal of Guidance, Control, and Dynamics , vol. 21,no. 2, pp. 193-207, March 1998.[17] M. Kelly, “An Introduction to Trajectory Optimization: How to Do Your Own Direct Collocation,”
SIAM Review , vol. 59, no. 4,pp. 849-904, 2017.[18] I.M. Ross, “A Primer on Pontryagin’s Principle in Optimal Control,” San Francisco, CA, USA: Collegiate Publishers, 2015.[19] J.T. Betts,
Practical Methods for Optimal Control and Estimation using Nonlinear Programming , Philadelphia, PA, USA:SIAM, 2010.[20] J. Nocedal, S.J. Wright,
Numerical Optimization , New York, NY, USA: Springer, 2006.[21] M. Szmuk, U. Eren, B. Açıkmeşe, “Successive Convexification for Mars 6-DoF Powered Descent Landing Guidance,” in
Proceedings of the AIAA SciTech Forum , Grapevine, TX, USA, January 2017.[22] M. Szmuk, B. Açıkmeşe, “Successive Convexification for 6-DoF Mars Rocket Powered Landing with Free-Final-Time,” in
Proceedings of the AIAA SciTech Forum , San Diego, CA, USA, January 2018.[23] D. Dueri, Y. Mao, Z. Mian, J. Ding, B. Açıkmeşe, “Trajectory Optimization with Inter-sample Obstacle Avoidance viaSuccessive Convexification,” in
Proceedings of the IEEE Conference on Decision and Control (CDC) , Melbourne, Australia,December 2017.[24] D. Malyuta, T.P. Reynolds, M. Szmuk, M. Mesbahi, B. Açıkmeşe, “Discretization Performance and Accuracy Analysis for thePowered Descent Guidance Problem,” in
Proceedings of the AIAA SciTech Forum , San Diego, CA, USA, January 2019.[25] Y. Mao, M. Szmuk, B. Açıkmeşe, “A Tutorial on Real-time Convex Optimization Based Guidance and Control for AerospaceApplications,” in
Proceedings of the IEEE American Control Conference , Milwaukee, WI, USA, June 2018.[26] Y. Mao, M. Szmuk, X. Xu, B. Açıkmeşe, “Successive Convexification: A Superlinearly Convergent Algorithm for Non-ConvexOptimal Control Problems,” arXiv:1804.06539 [math.OC] , February 2019.
27] M. Szmuk, T.P. Reynolds, B. Açıkmeşe, “Successive Convexification for Real-Time Six-Degree-of-Freedom Powered DescentGuidance with State-Triggered Constraints,”
Journal of Guidance, Control, and Dynamics , vol. 43, no. 8, pp. 1399-1413,August 2020.[28] T.P. Reynolds, M. Szmuk, D. Malyuta, M. Mesbahi, B. Açıkmeşe, “Dual Quaternion-Based Powered Descent Guidance withState-Triggered Constraints,”
Journal of Guidance, Control, and Dynamics , vol. 43, no. 9, pp. 1584-1599, September 2020.[29] M. Szmuk, D. Malyuta, T.P. Reynolds, M.S. Mceowen, B. Açıkmeşe, “Real-Time Quad-Rotor Path Planning Using ConvexOptimization and Compound State-Triggered Constraints,” arXiv:1902.09149 [math.OC] , February 2019.[30] M. Szmuk, T.P. Reynolds, B. Açıkmeşe, M. Mesbahi, “Successive Convexification for 6-DoF Powered Descent Guidance withCompound State-Triggered Constraints,” in
Proceedings of the AIAA SciTech Forum , San Diego, CA, USA, January 2019.[31] D. Malyuta, T.P. Reynolds, M. Szmuk, B. Ackimese, M. Mesbahi, “Fast Trajectory Optimization via Successive Convexificationfor Spacecraft Rendezvous with Integer Constraints,” arXiv:1906.04857 [math.OC] , June 2019.[32] M. Szmuk, “Successive Convexification & High Performance Feedback Control for Agile Flight,” PhD dissertation, Universityof Washington, 2019.[33] T.P. Reynolds, D. Malyuta, M. Mesbahi, B. Açıkmeşe, “A Real-Time Algorithm for Non-Convex Powered Descent Guidance,”in
Proceedings of the AIAA SciTech Forum , Orlando, FL, USA, January 2020.[34] P.C. Hughes,
Spacecraft Attitude Dynamics , Mineola, New York, USA: Dover Publications, 1986.[35] J.R. Wertz,
Spacecraft Attitude Determination and Control , Dordrecht, Holland: D. Reidel Publishing Company, 1978.[36] F.L. Markley, J.L. Crassidis,
Fundamentals of Spacecraft Attitude Determination and Control , New York, USA: Springer, 2014.[37] D. Hull, “Conversion of Optimal Control Problems into Parameter Optimization Problems,”
Journal of Guidance, Control, andDynamics , vol. 20, no. 1, pp. 57-60, January 1997.[38] Van Ryswyk, “Planet announces 50cm Skysat imagery, tasking dashboard and up to 12x Revisit,”
Planet Pulse , June 2020.[39] “NASA Earth Observatory: Earth Observing-1 (EO-1),” https://earthobservatory.nasa.gov/features/EO1/eo1_2.php , November 2000.[40] K. Shoemake, “Animating rotation with quaternion curves,” in
Proceedings of SIGGRAPH Computer Graphics , San Francisco,CA, USA, July 1985.[41] M. Grant, S. Boyd, “CVX: Matlab software for disciplined convex programming, version 2.0 beta,” http://cvxr.com/cvx ,September 2013.[42] M. Athan, P.L. Falb,
Optimal Control: An Introduction to the Theory and Its Applications , New York, USA: McGraw Hill, 1966.[43] A.E. Bryson JR, Y.C. Ho,
Applied Optimal Control , New York, NY, USA: Taylor & Francis, 1975.[44] B. Anderson, J. Moore, “Optimal Control: Linear Quadratic Methods”, Englewood Cliffs, NJ, USA: Prentice-Hall, 1990., New York, NY, USA: Taylor & Francis, 1975.[44] B. Anderson, J. Moore, “Optimal Control: Linear Quadratic Methods”, Englewood Cliffs, NJ, USA: Prentice-Hall, 1990.