A Convex Formulation of the Soft-Capture Problem

Ibrahima Sory Sow, Geordan Gutow, Howie Choset, Zachary Manchester *Ibrahima Sory Sow was supported by a Fellowship of the Belgian American Educational Foundation.Authors are with the Robotics Institute and Department of Mechanical Engineering, Carnegie Mellon University, 5000 Forbes Ave., Pittsburgh, PA 15213. {isow, ggutow, choset, zmanches}@andrew.cmu.edu
Abstract

We present a fast trajectory optimization algorithm for the soft capture of uncooperative tumbling space objects. Our algorithm generates safe, dynamically feasible, and minimum-fuel trajectories for a six-degree-of-freedom servicing spacecraft to achieve soft capture (near-zero relative velocity at contact) between predefined locations on the servicer spacecraft and target body. We solve a convex problem by enforcing a convex relaxation of the field-of-view constraint, followed by a sequential convex program correcting the trajectory for collision avoidance. The optimization problems can be solved with a standard second-order cone programming solver, making the algorithm both fast and practical for implementation in flight software. We demonstrate the performance and robustness of our algorithm in simulation over a range of object tumble rates up to 10/s.

I Introduction

Many envisioned autonomous on-orbit interactions, including active debris removal, refueling, and repair or upgrade missions, require an active chaser spacecraft to rendezvous and make contact with a passive or uncooperative target orbiting body [1]. A key requirement is soft capture: near-zero relative velocity at the time of contact. Rendezvous and docking operations have historically been limited to cooperative targets that communicate with the chaser and control their state to facilitate docking [2]. Such operations also require significant human interaction. However, capturing uncooperative or tumbling targets is also of significant interest. In such scenarios, human involvement is impractical due to the fast sense-plan-act loop required [3]. Therefore, the chaser needs to carry out agile maneuvers autonomously.

Refer to caption
Figure 1: The soft-capture maneuver exhibits four phases: 1) The chaser (beige) approaches the target (pink) while maintaining its line of sight. 2) Impulsive burns help avoid collision with the client (thrust vector in red). 3) Corrective burns align the capture end-effector with the target capture frame. 4) the chaser is within the capture ball and turns on its capture end-effector.

A soft-capture mission would start with a stand-off phase estimating the dynamic properties and tumble state of the target. Next, during the approach phase (0-50m), the chaser needs to achieve the necessary relative states to initiate the capture phase, when the capture interface is activated to grapple the target. A significant body of work has focused on using a free-floating robotic manipulator to grapple a slowly rotating target with minimal base disturbance [4]. Virgili-Llop and Romano [5] introduced a two-phase sequential-convex-programming (SCP) method for the capture of a tumbling target. For the approach phase, Sternberg and Miller [6] analyzed various parameterization schemes for creating fuel-efficient trajectories, whereas Stoneman and Lampariello [7] devised a computationally expensive approach for generating offline reference trajectories for online tracking. Albee et al. [8] experimented with offline-generated trajectory look-up tables for online parsing and robust tracking. Malyuta et al. [9] proposed an SCP-based strategy for six-degree-of-freedom docking maneuvers against a stationary target.

In this work, we focus on the approach phase, during which the majority of the mission ΔVΔ𝑉\Delta{V}roman_Δ italic_V is expended [3]. Viewed as a trajectory optimization problem, inherent nonconvexities in the constraints must be addressed. For instance, field of view constraints (direct line of sight with the target) have historically been relaxed through semidefinite programming (SDP) methods [10] [11] or through Convex-Concave programming approaches [12]. For collision-avoidance constraints, relaxations have relied on pre-determined keep-out zones or rotating hyperplanes requiring precise prior timing knowledge [13].

Our approach leverages recent advances in convex optimization, which have been applied to many complex aerospace guidance and control problems, such as the classical problem of powered soft landing [14]. Convex optimization offers deterministic guarantees given a desired solution accuracy [15] and has modest computational requirements, making it suitable for implementation onboard spacecraft.

We introduce a safe and efficient trajectory optimization algorithm for the chaser spacecraft to capture a target tumbling at up to 10/s{}^{\circ}/\mathrm{s}start_FLOATSUPERSCRIPT ∘ end_FLOATSUPERSCRIPT / roman_s. Our algorithm produces minimum-fuel trajectories that respect stringent and nonconvex safety constraints: target within field-of-view, collision avoidance, and near-zero relative velocity at the time of capture. We first solve a convex problem by enforcing a convex relaxation of the field-of-view constraints, followed by a sequential convex program correcting the trajectory for collision avoidance. The optimizations can be solved with a standard second-order cone programming solver, making the algorithm both fast and practical for implementation in flight software. We evaluate the convergence behavior and performance of the algorithm on 250 different target tumble rates (up to 10/s{}^{\circ}/\mathrm{s}start_FLOATSUPERSCRIPT ∘ end_FLOATSUPERSCRIPT / roman_s) and relative chaser orbits in simulation.

The paper proceeds as follows: Sec. II reviews relevant mathematics. In Sec. III the problem setup and dynamics are detailed. Sec. IV presents the nonconvex problem and Sec. V describes its convex relaxation. Sec. VI evaluates the computational, safety, and robustness performance of our algorithm, and Sec. VII summarizes our conclusions.

II Background

II-A Quaternion Algebra

We represent attitude with unit quaternions following the conventions in [16]. We define a quaternion q4:=[qsqvT]T𝑞superscript4assignsuperscriptdelimited-[]subscript𝑞𝑠superscriptsubscript𝑞𝑣𝑇𝑇q\in{\mathbb{R}}^{4}:=[q_{s}\;\;q_{v}^{T}]^{T}italic_q ∈ blackboard_R start_POSTSUPERSCRIPT 4 end_POSTSUPERSCRIPT := [ italic_q start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT italic_q start_POSTSUBSCRIPT italic_v end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_T end_POSTSUPERSCRIPT ] start_POSTSUPERSCRIPT italic_T end_POSTSUPERSCRIPT where qssubscript𝑞𝑠q_{s}\in{\mathbb{R}}italic_q start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT ∈ blackboard_R and qv3subscript𝑞𝑣superscript3q_{v}\in{\mathbb{R}}^{3}italic_q start_POSTSUBSCRIPT italic_v end_POSTSUBSCRIPT ∈ blackboard_R start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT are its scalar and vector parts, respectively. We define the matrix

L(q):=[qsqvTqvqsI+[qv]×]assign𝐿𝑞matrixsubscript𝑞𝑠superscriptsubscript𝑞𝑣𝑇subscript𝑞𝑣subscript𝑞𝑠𝐼superscriptdelimited-[]subscript𝑞𝑣\displaystyle L(q):=\begin{bmatrix}q_{s}\;\;&-q_{v}^{T}\\ q_{v}\;\;&q_{s}I+[q_{v}]^{\times}\end{bmatrix}italic_L ( italic_q ) := [ start_ARG start_ROW start_CELL italic_q start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT end_CELL start_CELL - italic_q start_POSTSUBSCRIPT italic_v end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_T end_POSTSUPERSCRIPT end_CELL end_ROW start_ROW start_CELL italic_q start_POSTSUBSCRIPT italic_v end_POSTSUBSCRIPT end_CELL start_CELL italic_q start_POSTSUBSCRIPT italic_s end_POSTSUBSCRIPT italic_I + [ italic_q start_POSTSUBSCRIPT italic_v end_POSTSUBSCRIPT ] start_POSTSUPERSCRIPT × end_POSTSUPERSCRIPT end_CELL end_ROW end_ARG ] (1)

to compute the quaternion product q2q1=L(q2)q1tensor-productsubscript𝑞2subscript𝑞1𝐿subscript𝑞2subscript𝑞1q_{2}\otimes q_{1}=L(q_{2})q_{1}italic_q start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT ⊗ italic_q start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT = italic_L ( italic_q start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT ) italic_q start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT. We use Q(q)𝑄𝑞Q(q)italic_Q ( italic_q ) to denote the rotation matrix equivalent to q𝑞qitalic_q.

II-B Differentiable Collision Metric

Collision-avoidance constraints yield a nonconvex feasible region that is time-varying and dependent on both orbit and attitude. Traditional collision-detection routines are inherently non-differentiable, making them hard to incorporate in optimization routines. In this work, we conservatively approximate the chaser and target geometries by their convex hulls through a set of linear inequalities and adopt DCOL [17], a differentiable collision formulation visualized in Fig. 2. DCOL solves a small convex optimization problem for the minimal inflation factor α𝛼\alphaitalic_α leading to an intersection between the two inflated convex hulls. The approach is differentiable by design, and collision is detected if α<1𝛼1\alpha<1italic_α < 1. DCOL allows fast computation of the Jacobian of the inflation factor α𝛼\alphaitalic_α with respect to the position and attitudes of both bodies [17], permitting its use for relaxations in sequential convex optimization.

Refer to caption
Figure 2: Left: Collision detection methods find a non-differentiable minimum distance between the 3D meshes. Right: DCOL [17] finds the minimum inflating factor α𝛼\alphaitalic_α between the convex hulls leading to an intersection with fast Jacobian computations.

III Relative Orbital Dynamics

Refer to caption
Figure 3: {O}𝑂\{O\}{ italic_O } is the reference frame. The target body frame {T}𝑇\{T\}{ italic_T }, rotating with an angular velocity vector ωtsubscript𝜔𝑡\omega_{t}italic_ω start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT, coincides with the center of {O}𝑂\{O\}{ italic_O }. Dtsubscript𝐷𝑡D_{t}italic_D start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT is the capture point on the target vehicle. The chaser body frame {C}𝐶\{C\}{ italic_C }, with relative position r𝑟ritalic_r and velocity v𝑣vitalic_v expressed in {O}𝑂\{O\}{ italic_O }, is centered at the center of mass of the chaser. u𝑢uitalic_u is the thrust vector. Dcsubscript𝐷𝑐D_{c}italic_D start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT is the chaser capture point.

We assume a passive (unactuated) and uncooperative target spacecraft in orbit around Earth. Its position coincides with the origin of a local-vertical-local-horizontal (LVLH) frame, or Hill frame, denoted {O}𝑂\{O\}{ italic_O }, where the X-axis points along the radial direction, the Z-axis along the orbital angular momentum vector, and the Y-axis completes the right-handed frame, pointing along the direction of motion. The attitude and angular rates of the target spacecraft make up the state xt(t)=[qtωt]subscript𝑥𝑡𝑡delimited-[]subscript𝑞𝑡subscript𝜔𝑡x_{t}(t)=[q_{t}\;\;\omega_{t}]italic_x start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT ( italic_t ) = [ italic_q start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT italic_ω start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT ] and are described through a frame {T}𝑇\{T\}{ italic_T } attached to its assumed center of mass as shown in Fig. 3. The tumbling trajectory follows Euler’s equations for rigid bodies and the quaternion kinematics equation:

xt=[qtωt],xt˙=[12L(qt)[0ωt]J1(ωt×Jωt)],formulae-sequencesubscript𝑥𝑡matrixsubscript𝑞𝑡subscript𝜔𝑡˙subscript𝑥𝑡matrix12𝐿subscript𝑞𝑡matrix0subscript𝜔𝑡superscript𝐽1subscript𝜔𝑡𝐽subscript𝜔𝑡x_{t}=\begin{bmatrix}q_{t}\\ \omega_{t}\end{bmatrix},\quad\dot{x_{t}}=\begin{bmatrix}\frac{1}{2}L(q_{t})% \begin{bmatrix}0\\ \omega_{t}\end{bmatrix}\\ J^{-1}(-\omega_{t}\times J\omega_{t})\end{bmatrix},italic_x start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT = [ start_ARG start_ROW start_CELL italic_q start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT end_CELL end_ROW start_ROW start_CELL italic_ω start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT end_CELL end_ROW end_ARG ] , over˙ start_ARG italic_x start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT end_ARG = [ start_ARG start_ROW start_CELL divide start_ARG 1 end_ARG start_ARG 2 end_ARG italic_L ( italic_q start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT ) [ start_ARG start_ROW start_CELL 0 end_CELL end_ROW start_ROW start_CELL italic_ω start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT end_CELL end_ROW end_ARG ] end_CELL end_ROW start_ROW start_CELL italic_J start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT ( - italic_ω start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT × italic_J italic_ω start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT ) end_CELL end_ROW end_ARG ] , (2)

where ωtsubscript𝜔𝑡\omega_{t}italic_ω start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT is the angular velocity of {T}𝑇\{T\}{ italic_T }, qtsubscript𝑞𝑡q_{t}italic_q start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT the rotation from {T}𝑇\{T\}{ italic_T } to {O}𝑂\{O\}{ italic_O }, and J𝐽Jitalic_J is the target’s inertia matrix.

We model the chaser spacecraft as a point mass with state x(t)=[r(t),v(t)]T=[rx(t),ry(t),rz(t),vx(t),vy(t),vz(t)]T𝑥𝑡superscript𝑟𝑡𝑣𝑡𝑇superscriptsubscript𝑟𝑥𝑡subscript𝑟𝑦𝑡subscript𝑟𝑧𝑡subscript𝑣𝑥𝑡subscript𝑣𝑦𝑡subscript𝑣𝑧𝑡𝑇x(t)=[r(t),v(t)]^{T}=[r_{x}(t),r_{y}(t),r_{z}(t),v_{x}(t),v_{y}(t),v_{z}(t)]^{T}italic_x ( italic_t ) = [ italic_r ( italic_t ) , italic_v ( italic_t ) ] start_POSTSUPERSCRIPT italic_T end_POSTSUPERSCRIPT = [ italic_r start_POSTSUBSCRIPT italic_x end_POSTSUBSCRIPT ( italic_t ) , italic_r start_POSTSUBSCRIPT italic_y end_POSTSUBSCRIPT ( italic_t ) , italic_r start_POSTSUBSCRIPT italic_z end_POSTSUBSCRIPT ( italic_t ) , italic_v start_POSTSUBSCRIPT italic_x end_POSTSUBSCRIPT ( italic_t ) , italic_v start_POSTSUBSCRIPT italic_y end_POSTSUBSCRIPT ( italic_t ) , italic_v start_POSTSUBSCRIPT italic_z end_POSTSUBSCRIPT ( italic_t ) ] start_POSTSUPERSCRIPT italic_T end_POSTSUPERSCRIPT and a thrust vector u(t)=[ux,uy,uz]T𝑢𝑡superscriptsubscript𝑢𝑥subscript𝑢𝑦subscript𝑢𝑧𝑇u(t)=[u_{x},u_{y},u_{z}]^{T}italic_u ( italic_t ) = [ italic_u start_POSTSUBSCRIPT italic_x end_POSTSUBSCRIPT , italic_u start_POSTSUBSCRIPT italic_y end_POSTSUBSCRIPT , italic_u start_POSTSUBSCRIPT italic_z end_POSTSUBSCRIPT ] start_POSTSUPERSCRIPT italic_T end_POSTSUPERSCRIPT. We decouple translation and attitude by assuming the presence of a low-level target-pointing attitude controller that can track any angular velocity reference ω(t)𝜔𝑡\omega(t)italic_ω ( italic_t ) with magnitude less than ωmaxsubscript𝜔𝑚𝑎𝑥\omega_{max}italic_ω start_POSTSUBSCRIPT italic_m italic_a italic_x end_POSTSUBSCRIPT. Given the short duration of the maneuver compared to the orbital period, the chaser thrust represents the most significant perturbation, and other orbital perturbations can be safely neglected. Under the assumption of a circular orbit for the target, the chaser motion can be modeled with the Hill-Clohessy-Wiltshire equations [18]

v˙xsubscript˙𝑣𝑥\displaystyle\dot{v}_{x}over˙ start_ARG italic_v end_ARG start_POSTSUBSCRIPT italic_x end_POSTSUBSCRIPT =3n2rx+2nvy+uxm,absent3superscript𝑛2subscript𝑟𝑥2𝑛subscript𝑣𝑦subscript𝑢𝑥𝑚\displaystyle=3n^{2}r_{x}+2nv_{y}+\frac{u_{x}}{m},= 3 italic_n start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_r start_POSTSUBSCRIPT italic_x end_POSTSUBSCRIPT + 2 italic_n italic_v start_POSTSUBSCRIPT italic_y end_POSTSUBSCRIPT + divide start_ARG italic_u start_POSTSUBSCRIPT italic_x end_POSTSUBSCRIPT end_ARG start_ARG italic_m end_ARG ,
v˙ysubscript˙𝑣𝑦\displaystyle\dot{v}_{y}over˙ start_ARG italic_v end_ARG start_POSTSUBSCRIPT italic_y end_POSTSUBSCRIPT =2nvx+uym,absent2𝑛subscript𝑣𝑥subscript𝑢𝑦𝑚\displaystyle=-2nv_{x}+\frac{u_{y}}{m},= - 2 italic_n italic_v start_POSTSUBSCRIPT italic_x end_POSTSUBSCRIPT + divide start_ARG italic_u start_POSTSUBSCRIPT italic_y end_POSTSUBSCRIPT end_ARG start_ARG italic_m end_ARG , (3)
v˙zsubscript˙𝑣𝑧\displaystyle\dot{v}_{z}over˙ start_ARG italic_v end_ARG start_POSTSUBSCRIPT italic_z end_POSTSUBSCRIPT =n2rz+uzm,absentsuperscript𝑛2subscript𝑟𝑧subscript𝑢𝑧𝑚\displaystyle=-n^{2}r_{z}+\frac{u_{z}}{m},= - italic_n start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_r start_POSTSUBSCRIPT italic_z end_POSTSUBSCRIPT + divide start_ARG italic_u start_POSTSUBSCRIPT italic_z end_POSTSUBSCRIPT end_ARG start_ARG italic_m end_ARG ,

which can be put in the standard matrix form

x˙(t)=Ax(t)+Bu(t)˙𝑥𝑡𝐴𝑥𝑡𝐵𝑢𝑡\dot{x}(t)=Ax(t)+Bu(t)over˙ start_ARG italic_x end_ARG ( italic_t ) = italic_A italic_x ( italic_t ) + italic_B italic_u ( italic_t ) (4)

where m𝑚mitalic_m is the mass of the chaser spacecraft, n=μ/a3𝑛𝜇superscript𝑎3n=\sqrt{\mu/a^{3}}italic_n = square-root start_ARG italic_μ / italic_a start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT end_ARG is the orbital mean motion of the target spacecraft, μ𝜇\muitalic_μ is the gravitational parameter for Earth, and a𝑎aitalic_a is the semi-major axis of the target orbit. The continuous dynamics equations are integrated to yield the discrete-time dynamics

xk+1=Adxk+Bduksubscript𝑥𝑘1subscript𝐴𝑑subscript𝑥𝑘subscript𝐵𝑑subscript𝑢𝑘x_{k+1}=A_{d}x_{k}+B_{d}u_{k}italic_x start_POSTSUBSCRIPT italic_k + 1 end_POSTSUBSCRIPT = italic_A start_POSTSUBSCRIPT italic_d end_POSTSUBSCRIPT italic_x start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT + italic_B start_POSTSUBSCRIPT italic_d end_POSTSUBSCRIPT italic_u start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT (5)

We assume that the chaser spacecraft initiates its proximity maneuvers from a passively stable relative orbit [19], i.e. a bounded relative orbit where the energies of the chaser and target orbits are matched. With a slightly different inclination and eccentricity from the target’s orbit, maintaining this orbit requires minimal ΔΔ\Deltaroman_ΔV. The set of initial conditions of the chaser is then restricted by

r˙y=2nrx,ry=2rx˙n,formulae-sequencesubscript˙𝑟𝑦2𝑛subscript𝑟𝑥subscript𝑟𝑦2˙subscript𝑟𝑥𝑛\dot{r}_{y}=-2nr_{x},\quad r_{y}=\frac{2\dot{r_{x}}}{n},over˙ start_ARG italic_r end_ARG start_POSTSUBSCRIPT italic_y end_POSTSUBSCRIPT = - 2 italic_n italic_r start_POSTSUBSCRIPT italic_x end_POSTSUBSCRIPT , italic_r start_POSTSUBSCRIPT italic_y end_POSTSUBSCRIPT = divide start_ARG 2 over˙ start_ARG italic_r start_POSTSUBSCRIPT italic_x end_POSTSUBSCRIPT end_ARG end_ARG start_ARG italic_n end_ARG , (6)

respectively canceling out the relative linear drift and offset of the relative orbit. Using the closed-form solution of the Clohessy-Wiltshire equations [18], the set of states lie on a parameterized “safe ellipsoid” with

A0=(rxn)2+rx2,B0=(rzn)2+rz2formulae-sequencesubscript𝐴0superscriptsubscript𝑟𝑥𝑛2superscriptsubscript𝑟𝑥2subscript𝐵0superscriptsubscript𝑟𝑧𝑛2superscriptsubscript𝑟𝑧2A_{0}=\sqrt{(\frac{r_{x}}{n})^{2}+r_{x}^{2}},\quad B_{0}=\sqrt{(\frac{r_{z}}{n% })^{2}+r_{z}^{2}}italic_A start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT = square-root start_ARG ( divide start_ARG italic_r start_POSTSUBSCRIPT italic_x end_POSTSUBSCRIPT end_ARG start_ARG italic_n end_ARG ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT + italic_r start_POSTSUBSCRIPT italic_x end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG , italic_B start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT = square-root start_ARG ( divide start_ARG italic_r start_POSTSUBSCRIPT italic_z end_POSTSUBSCRIPT end_ARG start_ARG italic_n end_ARG ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT + italic_r start_POSTSUBSCRIPT italic_z end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG (7)

where the radial, along-track, and cross-track dimensions of the ellipsoid are A0subscript𝐴0A_{0}italic_A start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT, 2A02subscript𝐴02A_{0}2 italic_A start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT, and B0subscript𝐵0B_{0}italic_B start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT, respectively. Combining (6) and (7) allows us to sample passively stable relative orbits.

The chaser predicts the future nominal motion of the target by integrating Eqs. (2) at a fixed timestep, where any state at t>0𝑡0t>0italic_t > 0 is accessed through cubic spline interpolation of the nominal trajectory. While we have chosen a spacecraft as a target, our method extends to any other passively tumbling space object.

IV The Soft-Capture Problem

The soft capture problem seeks state and control trajectories, x(t)𝑥𝑡x(t)italic_x ( italic_t ) and u(t)𝑢𝑡u(t)italic_u ( italic_t ), that drive the endpoint of the capture vector on the chaser spacecraft, Dcsubscript𝐷𝑐D_{c}italic_D start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT, to coincide with the endpoint of the capture vector on the target, Dtsubscript𝐷𝑡D_{t}italic_D start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT. The trajectories must minimize fuel consumption and satisfy mission-specific and safety constraints given the target state prediction xt(t)subscript𝑥𝑡𝑡x_{t}(t)italic_x start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT ( italic_t ). We pose this kinodynamic planning problem as a trajectory optimization problem of the form {mini} x_1:N, u_1:N-1ℓ_N (x_N) + ∑_k^N-1 ℓ_k (x_k, u_k) \addConstraintx_k+1 = A_d x_k + B_d u_k \addConstraintf(x_k, u_k)=0 \addConstraintg(x_k, u_k)≤0 where xksubscript𝑥𝑘x_{k}italic_x start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT and uksubscript𝑢𝑘u_{k}italic_u start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT are the decision variables, ksubscript𝑘\ell_{k}roman_ℓ start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT and Nsubscript𝑁\ell_{N}roman_ℓ start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT are the stage and terminal costs, and f𝑓fitalic_f and g𝑔gitalic_g specify equality and inequality constraints. If the objective and the constraints are convex in xk,uksubscript𝑥𝑘subscript𝑢𝑘x_{k},u_{k}italic_x start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT , italic_u start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT, the problem is a convex optimization problem, which can be solved to global optimality in polynomial time [15].

To begin with, we upper bound the thrust and velocity magnitudes to reflect thruster limitations and operational safety constraints:

uk2subscriptdelimited-∥∥subscript𝑢𝑘2\displaystyle\lVert u_{k}\rVert_{2}∥ italic_u start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT ∥ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT Umaxabsentsubscript𝑈𝑚𝑎𝑥\displaystyle\leq U_{max}≤ italic_U start_POSTSUBSCRIPT italic_m italic_a italic_x end_POSTSUBSCRIPT (8)
vk2\displaystyle\lVert{}v_{k}\rVert{}_{2}∥ italic_v start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT ∥ start_FLOATSUBSCRIPT 2 end_FLOATSUBSCRIPT vmaxabsentsubscript𝑣𝑚𝑎𝑥\displaystyle\leq v_{max}≤ italic_v start_POSTSUBSCRIPT italic_m italic_a italic_x end_POSTSUBSCRIPT (9)

Standard docking maneuvers require the chaser vehicle to approach the target in a cone-shaped corridor originating from the docking port of a relatively static target [1]. For our tumbling target scenario, we restrain the position of the chaser to be within a cone of half-angle θdocksubscript𝜃𝑑𝑜𝑐𝑘\theta_{dock}italic_θ start_POSTSUBSCRIPT italic_d italic_o italic_c italic_k end_POSTSUBSCRIPT (see Fig. 3) emanating from the capture point of the target spacecraft for the last Ndocksubscript𝑁𝑑𝑜𝑐𝑘N_{dock}italic_N start_POSTSUBSCRIPT italic_d italic_o italic_c italic_k end_POSTSUBSCRIPT timesteps. The resulting second-order cone constraint is

WXk2cTXktan(π2θdock)subscriptdelimited-∥∥𝑊subscript𝑋𝑘2superscript𝑐𝑇subscript𝑋𝑘𝜋2subscript𝜃𝑑𝑜𝑐𝑘\lVert WX_{k}\rVert_{2}\leq\frac{c^{T}X_{k}}{\tan(\frac{\pi}{2}-\theta_{dock})}∥ italic_W italic_X start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT ∥ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT ≤ divide start_ARG italic_c start_POSTSUPERSCRIPT italic_T end_POSTSUPERSCRIPT italic_X start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT end_ARG start_ARG roman_tan ( divide start_ARG italic_π end_ARG start_ARG 2 end_ARG - italic_θ start_POSTSUBSCRIPT italic_d italic_o italic_c italic_k end_POSTSUBSCRIPT ) end_ARG (10)

where Xksubscript𝑋𝑘X_{k}italic_X start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT is the position of the chaser relative to the target capture point and W𝑊Witalic_W and c𝑐citalic_c define the orientation and origin of the cone of the same target capture point

W𝑊\displaystyle Witalic_W =Rk[100010000],c=Rk[001]formulae-sequenceabsentsubscript𝑅𝑘matrix100010000𝑐subscript𝑅𝑘matrix001\displaystyle=R_{k}\begin{bmatrix}1&0&0\\ 0&1&0\\ 0&0&0\end{bmatrix},\quad c=R_{k}\begin{bmatrix}0\\ 0\\ 1\end{bmatrix}= italic_R start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT [ start_ARG start_ROW start_CELL 1 end_CELL start_CELL 0 end_CELL start_CELL 0 end_CELL end_ROW start_ROW start_CELL 0 end_CELL start_CELL 1 end_CELL start_CELL 0 end_CELL end_ROW start_ROW start_CELL 0 end_CELL start_CELL 0 end_CELL start_CELL 0 end_CELL end_ROW end_ARG ] , italic_c = italic_R start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT [ start_ARG start_ROW start_CELL 0 end_CELL end_ROW start_ROW start_CELL 0 end_CELL end_ROW start_ROW start_CELL 1 end_CELL end_ROW end_ARG ]

where Rk=Q(qt(kΔt))subscript𝑅𝑘𝑄subscript𝑞𝑡𝑘Δ𝑡R_{k}=Q(q_{t}(k\Delta t))italic_R start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT = italic_Q ( italic_q start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT ( italic_k roman_Δ italic_t ) ) is the attitude of the target at the timestep k𝑘kitalic_k.

The chaser’s attitude controller must maintain the line of sight of its optical sensors within some angular tolerance [1]. We ensure this by restricting the angle between line-of-sight vectors at successive timesteps to be less than ωmaxΔtsubscript𝜔𝑚𝑎𝑥Δ𝑡\omega_{max}\Delta titalic_ω start_POSTSUBSCRIPT italic_m italic_a italic_x end_POSTSUBSCRIPT roman_Δ italic_t:

rkTrk+1rkrk+1cos(ωmaxΔt)superscriptsubscript𝑟𝑘𝑇subscript𝑟𝑘1delimited-∥∥subscript𝑟𝑘delimited-∥∥subscript𝑟𝑘1subscript𝜔𝑚𝑎𝑥Δ𝑡\frac{r_{k}^{T}r_{k+1}}{\lVert r_{k}\rVert\lVert r_{k+1}\rVert}\geq\cos(\omega% _{max}\Delta t)divide start_ARG italic_r start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_T end_POSTSUPERSCRIPT italic_r start_POSTSUBSCRIPT italic_k + 1 end_POSTSUBSCRIPT end_ARG start_ARG ∥ italic_r start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT ∥ ∥ italic_r start_POSTSUBSCRIPT italic_k + 1 end_POSTSUBSCRIPT ∥ end_ARG ≥ roman_cos ( italic_ω start_POSTSUBSCRIPT italic_m italic_a italic_x end_POSTSUBSCRIPT roman_Δ italic_t ) (11)

Despite 0ωmaxΔtπ/20subscript𝜔𝑚𝑎𝑥Δ𝑡𝜋20\leq{}\omega_{max}\Delta t\leq{}\pi/20 ≤ italic_ω start_POSTSUBSCRIPT italic_m italic_a italic_x end_POSTSUBSCRIPT roman_Δ italic_t ≤ italic_π / 2 in our scenario, the constraint remains nonconvex. Note that this constraint, combined with Eq. (10), ensures correct alignment of the chaser for capture.

At the final time tfsubscript𝑡𝑓t_{f}italic_t start_POSTSUBSCRIPT italic_f end_POSTSUBSCRIPT, the endpoints of the capture vectors Dcsubscript𝐷𝑐D_{c}italic_D start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT and Dtsubscript𝐷𝑡D_{t}italic_D start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT must coincide (Eq. (12)) given the final predicted target’s attitude Qt(NΔt)subscript𝑄𝑡𝑁Δ𝑡Q_{t}(N\Delta t)italic_Q start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT ( italic_N roman_Δ italic_t ). The attitude controller ensures that Dcsubscript𝐷𝑐D_{c}italic_D start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT points at the target using the minimum-norm rotation increment at each timestep, so the attitude of the chaser is a position-dependent function, q(r)𝑞𝑟q(r)italic_q ( italic_r ), or its rotation matrix equivalent, Q(q(r))𝑄𝑞𝑟Q(q(r))italic_Q ( italic_q ( italic_r ) ). Thus, Q(q(rN))=Qt(NΔt)Qy(180)𝑄𝑞subscript𝑟𝑁subscript𝑄𝑡𝑁Δ𝑡subscript𝑄𝑦superscript180Q(q(r_{N}))=Q_{t}(N\Delta t)Q_{y}(-180^{\circ})italic_Q ( italic_q ( italic_r start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT ) ) = italic_Q start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT ( italic_N roman_Δ italic_t ) italic_Q start_POSTSUBSCRIPT italic_y end_POSTSUBSCRIPT ( - 180 start_POSTSUPERSCRIPT ∘ end_POSTSUPERSCRIPT ). To protect the integrity of the capture mechanism and the target, we bound the relative velocities between the target capture point and the chaser’s center of mass (Eq. (13))

rNQt(NΔt)(Dc+Dt)2subscriptdelimited-∥∥subscript𝑟𝑁subscript𝑄𝑡𝑁Δ𝑡subscript𝐷𝑐subscript𝐷𝑡2\displaystyle\lVert r_{N}-Q_{t}(N\Delta t)(D_{c}+D_{t})\rVert_{2}∥ italic_r start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT - italic_Q start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT ( italic_N roman_Δ italic_t ) ( italic_D start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT + italic_D start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT ) ∥ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT ϵpabsentsubscriptitalic-ϵ𝑝\displaystyle\leq\epsilon_{p}≤ italic_ϵ start_POSTSUBSCRIPT italic_p end_POSTSUBSCRIPT (12)
vNωt((N1)Δt)×Dt2subscriptdelimited-∥∥subscript𝑣𝑁subscript𝜔𝑡𝑁1Δ𝑡subscript𝐷𝑡2\displaystyle\lVert v_{N}-\omega_{t}((N-1)\Delta t)\times D_{t}\rVert_{2}∥ italic_v start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT - italic_ω start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT ( ( italic_N - 1 ) roman_Δ italic_t ) × italic_D start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT ∥ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT ϵvabsentsubscriptitalic-ϵ𝑣\displaystyle\leq\epsilon_{v}≤ italic_ϵ start_POSTSUBSCRIPT italic_v end_POSTSUBSCRIPT (13)

with ϵpsubscriptitalic-ϵ𝑝\epsilon_{p}italic_ϵ start_POSTSUBSCRIPT italic_p end_POSTSUBSCRIPT and ϵvsubscriptitalic-ϵ𝑣\epsilon_{v}italic_ϵ start_POSTSUBSCRIPT italic_v end_POSTSUBSCRIPT respectively the terminal position and velocity tolerances, guaranteeing the inequalities achievement while softening the constraints for the optimization solver.

The chaser must not collide with the target spacecraft during the trajectory. We use the collision metric described in Sec. II-B and require k:αk>αmin>1:for-all𝑘subscript𝛼𝑘subscript𝛼𝑚𝑖𝑛1\forall k:\alpha_{k}>\alpha_{min}>1∀ italic_k : italic_α start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT > italic_α start_POSTSUBSCRIPT italic_m italic_i italic_n end_POSTSUBSCRIPT > 1, ensuring a minimal separation distance between the convex hulls. The collision-avoidance constraint at timestep k𝑘kitalic_k is

α(rk,qk(rk),qt(kΔt))=α(rk,qt(kΔt))>αmin𝛼subscript𝑟𝑘subscript𝑞𝑘subscript𝑟𝑘subscript𝑞𝑡𝑘Δ𝑡𝛼subscript𝑟𝑘subscript𝑞𝑡𝑘Δ𝑡subscript𝛼𝑚𝑖𝑛\alpha(r_{k},q_{k}(r_{k}),q_{t}(k\Delta t))=\alpha(r_{k},q_{t}(k\Delta t))>% \alpha_{min}italic_α ( italic_r start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT , italic_q start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT ( italic_r start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT ) , italic_q start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT ( italic_k roman_Δ italic_t ) ) = italic_α ( italic_r start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT , italic_q start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT ( italic_k roman_Δ italic_t ) ) > italic_α start_POSTSUBSCRIPT italic_m italic_i italic_n end_POSTSUBSCRIPT (14)

To reduce fuel consumption, we minimize the sum of the L1 norms of the thrust vectors. The L1 norm of the control is proportional to fuel consumption and a better fuel minimizer than the (perhaps more common) L2 norm [20]. In addition, the reaction control system of the chaser employs impulsive burns, where the thrust magnitude is close either to zero or to the maximum limit [21]. The L1 norm encourages such “bang-off-bang” solutions. With this final addition, we fully define Problem 1: {mdframed} Problem 1: Nonconvex Soft Capture Problem {mini} x_k, u_k∑_k^N-1 ∥ u_k ∥_1 \addConstraint(5), (8), (9), (10) \addConstraint(11), (12), (13), (14)

Problem 1 is a nonconvex trajectory optimization problem where nonconvexities arise from the field-of-view and collision-avoidance constraints. Nonlinear programming solvers, such as IPOPT [22] and SNOPT [23], could potentially solve Problem 1 directly. However, such solvers lack time-complexity guarantees and converge to the nearest local minimum to some initial guess, complicating their real-time implementation.

Refer to caption
Figure 4: Illustration of the final algorithm. The future trajectory of the target qt(t),ωt(t)subscript𝑞𝑡𝑡subscript𝜔𝑡𝑡q_{t}(t),\omega_{t}(t)italic_q start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT ( italic_t ) , italic_ω start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT ( italic_t ) is predicted using the current estimate of the state of the target qt0,ωt0superscriptsubscript𝑞𝑡0superscriptsubscript𝜔𝑡0q_{t}^{0},\omega_{t}^{0}italic_q start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 0 end_POSTSUPERSCRIPT , italic_ω start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 0 end_POSTSUPERSCRIPT. The prediction, along the current state of the chaser r0,v0superscript𝑟0superscript𝑣0r^{0},v^{0}italic_r start_POSTSUPERSCRIPT 0 end_POSTSUPERSCRIPT , italic_v start_POSTSUPERSCRIPT 0 end_POSTSUPERSCRIPT, is input to an initial convex optimization problem without collision-avoidance constraints to obtain the first iterate of the trajectory x1,u1superscript𝑥1superscript𝑢1x^{1},u^{1}italic_x start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT , italic_u start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT for the sequential convex programming (SCP) pipeline. The SCP linearizes the collision constraints and iteratively solves convex subproblems until the safety criteria k:αk>1:for-all𝑘subscript𝛼𝑘1\forall k:\alpha_{k}>1∀ italic_k : italic_α start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT > 1 is met to yield xsol,usolsuperscript𝑥𝑠𝑜𝑙superscript𝑢𝑠𝑜𝑙x^{sol},u^{sol}italic_x start_POSTSUPERSCRIPT italic_s italic_o italic_l end_POSTSUPERSCRIPT , italic_u start_POSTSUPERSCRIPT italic_s italic_o italic_l end_POSTSUPERSCRIPT. If the SCP reaches a set maximum number of iterations, the algorithm reports infeasibility.

V Convex Trajectory Optimization

In this section, we convexify Problem 1 by relaxing the field-of-view (11) and collision-avoidance (14) constraints.

V-A Field Of View

We start by introducing an artificial bound constraint on the position vector rmaxsubscript𝑟𝑚𝑎𝑥r_{max}italic_r start_POSTSUBSCRIPT italic_m italic_a italic_x end_POSTSUBSCRIPT and lifting the dimensions of the state vector x𝑥xitalic_x by introducing the slack variable ρksubscript𝜌𝑘\rho_{k}italic_ρ start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT

0ρkrmax,rk2ρk.formulae-sequence0subscript𝜌𝑘subscript𝑟𝑚𝑎𝑥subscriptnormsubscript𝑟𝑘2subscript𝜌𝑘0\leq\rho_{k}\leq r_{max},\quad\left\|r_{k}\right\|_{2}\leq\rho_{k}.0 ≤ italic_ρ start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT ≤ italic_r start_POSTSUBSCRIPT italic_m italic_a italic_x end_POSTSUBSCRIPT , ∥ italic_r start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT ∥ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT ≤ italic_ρ start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT . (15)

Replacing the norms in (11) by ρ𝜌\rhoitalic_ρ and multiplying by -1 gives

rkTrk+1ρkρk+1cos(ωmaxΔt).superscriptsubscript𝑟𝑘𝑇subscript𝑟𝑘1subscript𝜌𝑘subscript𝜌𝑘1subscript𝜔𝑚𝑎𝑥Δ𝑡r_{k}^{T}r_{k+1}\leq-\rho_{k}\rho_{k+1}\cos(\omega_{max}\Delta t).italic_r start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_T end_POSTSUPERSCRIPT italic_r start_POSTSUBSCRIPT italic_k + 1 end_POSTSUBSCRIPT ≤ - italic_ρ start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT italic_ρ start_POSTSUBSCRIPT italic_k + 1 end_POSTSUBSCRIPT roman_cos ( italic_ω start_POSTSUBSCRIPT italic_m italic_a italic_x end_POSTSUBSCRIPT roman_Δ italic_t ) . (16)

Minimizing ρksubscript𝜌𝑘\rho_{k}italic_ρ start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT in the objective guarantees a tight relaxation [14]. We tackle the bilinearity by first noting that

rk2ρksubscriptnormsubscript𝑟𝑘2subscript𝜌𝑘\displaystyle\left\|r_{k}\right\|_{2}\leq\rho_{k}∥ italic_r start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT ∥ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT ≤ italic_ρ start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT rkTrkρk2,absentsuperscriptsubscript𝑟𝑘𝑇subscript𝑟𝑘superscriptsubscript𝜌𝑘2\displaystyle\implies{}r_{k}^{T}r_{k}\leq\rho_{k}^{2},⟹ italic_r start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_T end_POSTSUPERSCRIPT italic_r start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT ≤ italic_ρ start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT , (17)
rk+12ρk+1subscriptnormsubscript𝑟𝑘12subscript𝜌𝑘1\displaystyle\left\|r_{k+1}\right\|_{2}\leq\rho_{k+1}∥ italic_r start_POSTSUBSCRIPT italic_k + 1 end_POSTSUBSCRIPT ∥ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT ≤ italic_ρ start_POSTSUBSCRIPT italic_k + 1 end_POSTSUBSCRIPT rk+1Trk+1ρk+12.absentsuperscriptsubscript𝑟𝑘1𝑇subscript𝑟𝑘1superscriptsubscript𝜌𝑘12\displaystyle\implies{}r_{k+1}^{T}r_{k+1}\leq\rho_{k+1}^{2}.⟹ italic_r start_POSTSUBSCRIPT italic_k + 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_T end_POSTSUPERSCRIPT italic_r start_POSTSUBSCRIPT italic_k + 1 end_POSTSUBSCRIPT ≤ italic_ρ start_POSTSUBSCRIPT italic_k + 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT . (18)

Adding the resulting terms (Eqs. (17), (18)) respectively in both sides of Eq. (16) allows the expansion of a quadratic form

rkTrk+rkTrk+1+rkTrk+1ρk2ρkρk+1cos(ϕ)+ρk+12,superscriptsubscript𝑟𝑘𝑇subscript𝑟𝑘superscriptsubscript𝑟𝑘𝑇subscript𝑟𝑘1superscriptsubscript𝑟𝑘𝑇subscript𝑟𝑘1superscriptsubscript𝜌𝑘2subscript𝜌𝑘subscript𝜌𝑘1italic-ϕsuperscriptsubscript𝜌𝑘12r_{k}^{T}r_{k}+r_{k}^{T}r_{k+1}+r_{k}^{T}r_{k+1}\leq\rho_{k}^{2}-\rho_{k}\rho_% {k+1}\cos(\phi)+\rho_{k+1}^{2},italic_r start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_T end_POSTSUPERSCRIPT italic_r start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT + italic_r start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_T end_POSTSUPERSCRIPT italic_r start_POSTSUBSCRIPT italic_k + 1 end_POSTSUBSCRIPT + italic_r start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_T end_POSTSUPERSCRIPT italic_r start_POSTSUBSCRIPT italic_k + 1 end_POSTSUBSCRIPT ≤ italic_ρ start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT - italic_ρ start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT italic_ρ start_POSTSUBSCRIPT italic_k + 1 end_POSTSUBSCRIPT roman_cos ( italic_ϕ ) + italic_ρ start_POSTSUBSCRIPT italic_k + 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ,

where ϕ=ωmaxΔt>0italic-ϕsubscript𝜔𝑚𝑎𝑥Δ𝑡0\phi=\omega_{max}\Delta t>0italic_ϕ = italic_ω start_POSTSUBSCRIPT italic_m italic_a italic_x end_POSTSUBSCRIPT roman_Δ italic_t > 0. This can be rewritten in matrix form

zkTHzkykTSyk,superscriptsubscript𝑧𝑘𝑇𝐻subscript𝑧𝑘superscriptsubscript𝑦𝑘𝑇𝑆subscript𝑦𝑘z_{k}^{T}Hz_{k}\leq y_{k}^{T}Sy_{k},italic_z start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_T end_POSTSUPERSCRIPT italic_H italic_z start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT ≤ italic_y start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_T end_POSTSUPERSCRIPT italic_S italic_y start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT , (19)

with

zksubscript𝑧𝑘\displaystyle z_{k}italic_z start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT =[rkrk+1],H=[I3x30.5I3x30.5I3x3I3x3]>0,formulae-sequenceabsentmatrixsubscript𝑟𝑘subscript𝑟𝑘1𝐻matrixsubscript𝐼3𝑥30.5subscript𝐼3𝑥30.5subscript𝐼3𝑥3subscript𝐼3𝑥30\displaystyle=\begin{bmatrix}r_{k}\\ r_{k+1}\end{bmatrix},\quad H=\begin{bmatrix}I_{3x3}&0.5I_{3x3}\\ 0.5I_{3x3}&I_{3x3}\end{bmatrix}>0,= [ start_ARG start_ROW start_CELL italic_r start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT end_CELL end_ROW start_ROW start_CELL italic_r start_POSTSUBSCRIPT italic_k + 1 end_POSTSUBSCRIPT end_CELL end_ROW end_ARG ] , italic_H = [ start_ARG start_ROW start_CELL italic_I start_POSTSUBSCRIPT 3 italic_x 3 end_POSTSUBSCRIPT end_CELL start_CELL 0.5 italic_I start_POSTSUBSCRIPT 3 italic_x 3 end_POSTSUBSCRIPT end_CELL end_ROW start_ROW start_CELL 0.5 italic_I start_POSTSUBSCRIPT 3 italic_x 3 end_POSTSUBSCRIPT end_CELL start_CELL italic_I start_POSTSUBSCRIPT 3 italic_x 3 end_POSTSUBSCRIPT end_CELL end_ROW end_ARG ] > 0 ,
yksubscript𝑦𝑘\displaystyle y_{k}italic_y start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT =[ρkρk+1],S=[10.5cos(ϕ)0.5cos(ϕ)1]>0.formulae-sequenceabsentmatrixsubscript𝜌𝑘subscript𝜌𝑘1𝑆matrix10.5italic-ϕ0.5italic-ϕ10\displaystyle=\begin{bmatrix}\rho_{k}\\ \rho_{k+1}\end{bmatrix},\quad S=\begin{bmatrix}1&-0.5\cos(\phi)\\ -0.5\cos(\phi)&1\end{bmatrix}>0.= [ start_ARG start_ROW start_CELL italic_ρ start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT end_CELL end_ROW start_ROW start_CELL italic_ρ start_POSTSUBSCRIPT italic_k + 1 end_POSTSUBSCRIPT end_CELL end_ROW end_ARG ] , italic_S = [ start_ARG start_ROW start_CELL 1 end_CELL start_CELL - 0.5 roman_cos ( italic_ϕ ) end_CELL end_ROW start_ROW start_CELL - 0.5 roman_cos ( italic_ϕ ) end_CELL start_CELL 1 end_CELL end_ROW end_ARG ] > 0 .

As H𝐻Hitalic_H and S𝑆Sitalic_S are both symmetric and positive definite (since ϕ>0italic-ϕ0\phi>0italic_ϕ > 0), we can simplify Eq. (19)

Hzk2Syk2superscriptnorm𝐻subscript𝑧𝑘2superscriptnorm𝑆subscript𝑦𝑘2\left\|\sqrt{H}z_{k}\right\|^{2}\leq\left\|\sqrt{S}y_{k}\right\|^{2}∥ square-root start_ARG italic_H end_ARG italic_z start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT ∥ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ≤ ∥ square-root start_ARG italic_S end_ARG italic_y start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT ∥ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT (20)

and by taking out the square

HzkSyknorm𝐻subscript𝑧𝑘norm𝑆subscript𝑦𝑘\left\|\sqrt{H}z_{k}\right\|\leq\left\|\sqrt{S}y_{k}\right\|∥ square-root start_ARG italic_H end_ARG italic_z start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT ∥ ≤ ∥ square-root start_ARG italic_S end_ARG italic_y start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT ∥ (21)

The eigenvalues of S𝑆Sitalic_S can be computed in closed-form and its smallest one is λmin(S)=1cos(ϕ)subscript𝜆𝑚𝑖𝑛𝑆1𝑐𝑜𝑠italic-ϕ\lambda_{min}(S)=1-cos(\phi)italic_λ start_POSTSUBSCRIPT italic_m italic_i italic_n end_POSTSUBSCRIPT ( italic_S ) = 1 - italic_c italic_o italic_s ( italic_ϕ ). Subsequently, λmin(S)=1cos(ϕ)subscript𝜆𝑚𝑖𝑛𝑆1𝑐𝑜𝑠italic-ϕ\lambda_{min}(\sqrt{S})=\sqrt{1-cos(\phi)}italic_λ start_POSTSUBSCRIPT italic_m italic_i italic_n end_POSTSUBSCRIPT ( square-root start_ARG italic_S end_ARG ) = square-root start_ARG 1 - italic_c italic_o italic_s ( italic_ϕ ) end_ARG and we can write

Hzkλmin(S)yk2=1cos(ϕ)yk2.norm𝐻subscript𝑧𝑘subscript𝜆𝑚𝑖𝑛𝑆subscriptnormsubscript𝑦𝑘21𝑐𝑜𝑠italic-ϕsubscriptnormsubscript𝑦𝑘2\left\|\sqrt{H}z_{k}\right\|\leq\lambda_{min}(\sqrt{S})\left\|y_{k}\right\|_{2% }=\sqrt{1-cos(\phi)}\left\|y_{k}\right\|_{2}.∥ square-root start_ARG italic_H end_ARG italic_z start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT ∥ ≤ italic_λ start_POSTSUBSCRIPT italic_m italic_i italic_n end_POSTSUBSCRIPT ( square-root start_ARG italic_S end_ARG ) ∥ italic_y start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT ∥ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT = square-root start_ARG 1 - italic_c italic_o italic_s ( italic_ϕ ) end_ARG ∥ italic_y start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT ∥ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT . (22)

Consecutive position vectors in time slightly differ in magnitude, ρkρk+1subscript𝜌𝑘subscript𝜌𝑘1\rho_{k}\approx\rho_{k+1}italic_ρ start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT ≈ italic_ρ start_POSTSUBSCRIPT italic_k + 1 end_POSTSUBSCRIPT, and can be used to approximate yk2subscriptnormsubscript𝑦𝑘2\left\|y_{k}\right\|_{2}∥ italic_y start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT ∥ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT

yk2=ρk2+ρk+1212(ρk+ρk+1),subscriptnormsubscript𝑦𝑘2superscriptsubscript𝜌𝑘2superscriptsubscript𝜌𝑘1212subscript𝜌𝑘subscript𝜌𝑘1\left\|y_{k}\right\|_{2}=\sqrt{\rho_{k}^{2}+\rho_{k+1}^{2}}\approx\frac{1}{% \sqrt{2}}(\rho_{k}+\rho_{k+1}),∥ italic_y start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT ∥ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT = square-root start_ARG italic_ρ start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT + italic_ρ start_POSTSUBSCRIPT italic_k + 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG ≈ divide start_ARG 1 end_ARG start_ARG square-root start_ARG 2 end_ARG end_ARG ( italic_ρ start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT + italic_ρ start_POSTSUBSCRIPT italic_k + 1 end_POSTSUBSCRIPT ) ,

which leads to

Hzk1cos(ϕ)2(ρk+ρk+1).norm𝐻subscript𝑧𝑘1𝑐𝑜𝑠italic-ϕ2subscript𝜌𝑘subscript𝜌𝑘1\left\|\sqrt{H}z_{k}\right\|\leq\sqrt{\frac{1-cos(\phi)}{2}}(\rho_{k}+\rho_{k+% 1}).∥ square-root start_ARG italic_H end_ARG italic_z start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT ∥ ≤ square-root start_ARG divide start_ARG 1 - italic_c italic_o italic_s ( italic_ϕ ) end_ARG start_ARG 2 end_ARG end_ARG ( italic_ρ start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT + italic_ρ start_POSTSUBSCRIPT italic_k + 1 end_POSTSUBSCRIPT ) .

Finally, with the trigonometric relation sin2(ϕ)=1cos(ϕ)2𝑠𝑖superscript𝑛2italic-ϕ1𝑐𝑜𝑠italic-ϕ2sin^{2}(\phi)=\frac{1-cos(\phi)}{2}italic_s italic_i italic_n start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ( italic_ϕ ) = divide start_ARG 1 - italic_c italic_o italic_s ( italic_ϕ ) end_ARG start_ARG 2 end_ARG and 1cos(ϕ)>01𝑐𝑜𝑠italic-ϕ01-cos(\phi)>01 - italic_c italic_o italic_s ( italic_ϕ ) > 0, we derive the relaxed convex second-order cone field-of-view constraints:

Hzksin(ϕ2)(ρk+ρk+1).norm𝐻subscript𝑧𝑘𝑠𝑖𝑛italic-ϕ2subscript𝜌𝑘subscript𝜌𝑘1\left\|\sqrt{H}z_{k}\right\|\leq sin(\frac{\phi}{2})(\rho_{k}+\rho_{k+1}).∥ square-root start_ARG italic_H end_ARG italic_z start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT ∥ ≤ italic_s italic_i italic_n ( divide start_ARG italic_ϕ end_ARG start_ARG 2 end_ARG ) ( italic_ρ start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT + italic_ρ start_POSTSUBSCRIPT italic_k + 1 end_POSTSUBSCRIPT ) . (23)

V-B Collision Avoidance

We leverage the fast computation of Jacobians of the collision detection factor αksubscript𝛼𝑘\alpha_{k}italic_α start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT (Sec. II-B) with respect to rksubscript𝑟𝑘r_{k}italic_r start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT to linearize the nonconvex collision-avoidance constraints about a reference rkisubscriptsuperscript𝑟𝑖𝑘r^{i}_{k}italic_r start_POSTSUPERSCRIPT italic_i end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT, with the superscript i𝑖iitalic_i denoting the iteration number. We further add a positive slack variable sksubscript𝑠𝑘s_{k}italic_s start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT, relaxing the inequality and ensuring its feasibility:

αk(rki+δrk)subscript𝛼𝑘subscriptsuperscript𝑟𝑖𝑘𝛿subscript𝑟𝑘\displaystyle\alpha_{k}(r^{i}_{k}+\delta r_{k})italic_α start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT ( italic_r start_POSTSUPERSCRIPT italic_i end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT + italic_δ italic_r start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT ) α(rki,qt(kΔt))+Jkiδrk+sk>αmin,absent𝛼subscriptsuperscript𝑟𝑖𝑘subscript𝑞𝑡𝑘Δ𝑡superscriptsubscript𝐽𝑘𝑖𝛿subscript𝑟𝑘subscript𝑠𝑘subscript𝛼𝑚𝑖𝑛\displaystyle\approx\alpha(r^{i}_{k},q_{t}(k\Delta t))+J_{k}^{i}\delta r_{k}+s% _{k}>\alpha_{min},≈ italic_α ( italic_r start_POSTSUPERSCRIPT italic_i end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT , italic_q start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT ( italic_k roman_Δ italic_t ) ) + italic_J start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_i end_POSTSUPERSCRIPT italic_δ italic_r start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT + italic_s start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT > italic_α start_POSTSUBSCRIPT italic_m italic_i italic_n end_POSTSUBSCRIPT , (24)

where

Jki=αr|rki+αq(r)|rkiq(r)r|rki,superscriptsubscript𝐽𝑘𝑖evaluated-at𝛼𝑟subscriptsuperscript𝑟𝑖𝑘evaluated-atevaluated-at𝛼𝑞𝑟subscriptsuperscript𝑟𝑖𝑘𝑞𝑟𝑟subscriptsuperscript𝑟𝑖𝑘J_{k}^{i}=\frac{\partial\alpha}{\partial r}\bigg{|}_{r^{i}_{k}}+\frac{\partial% \alpha}{\partial q(r)}\bigg{|}_{r^{i}_{k}}\frac{\partial q(r)}{\partial r}% \bigg{|}_{r^{i}_{k}},italic_J start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_i end_POSTSUPERSCRIPT = divide start_ARG ∂ italic_α end_ARG start_ARG ∂ italic_r end_ARG | start_POSTSUBSCRIPT italic_r start_POSTSUPERSCRIPT italic_i end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT end_POSTSUBSCRIPT + divide start_ARG ∂ italic_α end_ARG start_ARG ∂ italic_q ( italic_r ) end_ARG | start_POSTSUBSCRIPT italic_r start_POSTSUPERSCRIPT italic_i end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT end_POSTSUBSCRIPT divide start_ARG ∂ italic_q ( italic_r ) end_ARG start_ARG ∂ italic_r end_ARG | start_POSTSUBSCRIPT italic_r start_POSTSUPERSCRIPT italic_i end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT end_POSTSUBSCRIPT ,

We penalize sksubscript𝑠𝑘s_{k}italic_s start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT in the objective function with a penalty weight ψ𝜓\psiitalic_ψ as a measure of constraint infeasibility. For sufficiently large values of ψ𝜓\psiitalic_ψ, the objective term ψsk𝜓subscript𝑠𝑘\psi s_{k}italic_ψ italic_s start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT defines an exact penalty function [24] and can be minimized by solving a sequence of subproblems.

V-C Final Algorithm

The final algorithm is depicted in Fig. 4. With the chaser state and target trajectory prediction as input, we solve an initial convex second-order-cone program (SOCP) for x,u𝑥𝑢x,uitalic_x , italic_u with the convexified field of view constraint but without the collision-avoidance constraints (Problem 2). The near-feasible and near-optimal trajectory is then used to initialize a sequential convex subroutine solving for the corrections δx,δu𝛿𝑥𝛿𝑢\delta x,\delta uitalic_δ italic_x , italic_δ italic_u where the collision-avoidance constraint is added back and successively linearized along the previous solution (Eq. (24)), yielding Problem 3. We then update the current solution i𝑖iitalic_i

xki+1=xki+δxk,uki+1=uki+δukformulae-sequencesuperscriptsubscript𝑥𝑘𝑖1superscriptsubscript𝑥𝑘𝑖𝛿subscript𝑥𝑘superscriptsubscript𝑢𝑘𝑖1superscriptsubscript𝑢𝑘𝑖𝛿subscript𝑢𝑘\displaystyle x_{k}^{i+1}=x_{k}^{i}+\delta x_{k},\quad u_{k}^{i+1}=u_{k}^{i}+% \delta u_{k}italic_x start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_i + 1 end_POSTSUPERSCRIPT = italic_x start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_i end_POSTSUPERSCRIPT + italic_δ italic_x start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT , italic_u start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_i + 1 end_POSTSUPERSCRIPT = italic_u start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_i end_POSTSUPERSCRIPT + italic_δ italic_u start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT (25)

This iterative convexification and correction scheme is repeated until a strictly safe trajectory k:αk>1:for-all𝑘subscript𝛼𝑘1\forall k:\alpha_{k}>1∀ italic_k : italic_α start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT > 1 has been achieved as illustrated in Fig. 4. This choice of convergence criteria allows successful trajectories to go below the soft conservative limit αminsubscript𝛼𝑚𝑖𝑛\alpha_{min}italic_α start_POSTSUBSCRIPT italic_m italic_i italic_n end_POSTSUBSCRIPT.

The slack variable sksubscript𝑠𝑘s_{k}italic_s start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT used in the collision-avoidance constraints ensures the feasibility of the subproblems and the exact penalty term in the objective function will succcessively eliminate any residual infeasibility. We report failure if Problem 2 is infeasible or if the SCP remains infeasible after the maximum allowable number of iterations (Problem 3), i.e. the solver has converged to a non-zero constraint violation.

While not guaranteed to converge to a feasible solution, the algorithm is an efficient heuristic and can quickly detect an infeasible problem. The collision-avoidance constraints are applied in an uncrowded space where only the chaser and target are at risk of collision. Thus, as we show in Section VI, only small corrections are required from Problem 3 when the solution to Problem 2 is in collision. This permits the omission of a trust region.

{mdframed}

Problem 2: Convex Soft Capture Problem without collision-avoidance constraints {mini} x_k, u_k, ρ_k∑_k^N-1 ∥ u_k ∥_1 + ∑_k^N ρ_k \addConstraint(5), (8), (9), (10),(12) \addConstraint(13), (15), (23)

{mdframed}

Problem 3: Convexified Soft Capture Sub-Problem with collision-avoidance constraints {mini} δx_k, δu_k, ρ_k, s_kγ∑_k^N-1 ∥ u^i + δu_k ∥_1 + ∑_k^N (ρ_k + ψs_k) \addConstraint(5), (8), (9), (10),(12) \addConstraint(13), (15), (23), (24)

Refer to caption
Figure 5: The initial trajectory (left) resulting from solving Problem 2 is infeasible. The SCP procedure (Problem 3) adjusts the “tail” of the trajectory to avoid collision with the target (right).

VI Numerical experiments

We evaluate our algorithm through two sets of experiments: an illustration of the corrective behavior of the SCP to satisfy the collision-avoidance constraint and an analysis of the computational and robustness properties of the algorithm for 250 pairs of sampled target initial tumble rates and chaser relative orbits. Numerical experiments were conducted on a laptop equipped with an Intel i7-1165G7 2.80GHz CPU and 16 GB of RAM using the optimization solver MOSEK [25], whose solve times are reported. All experiments were conducted using the Julia programming language. Parameters used for experiments are listed in Table I. The capture mechanism of the chaser is abstracted through the terminal position tolerance ϵpsubscriptitalic-ϵ𝑝\epsilon_{p}italic_ϵ start_POSTSUBSCRIPT italic_p end_POSTSUBSCRIPT. The target’s inertia matrix is:

J=[5.890560.00.00.011.44620.2335160.00.23351611.5365]𝐽matrix5.890560.00.00.011.44620.2335160.00.23351611.5365J=\begin{bmatrix}5.89056&0.0&0.0\\ 0.0&11.4462&0.233516\\ 0.0&0.233516&11.5365\\ \end{bmatrix}italic_J = [ start_ARG start_ROW start_CELL 5.89056 end_CELL start_CELL 0.0 end_CELL start_CELL 0.0 end_CELL end_ROW start_ROW start_CELL 0.0 end_CELL start_CELL 11.4462 end_CELL start_CELL 0.233516 end_CELL end_ROW start_ROW start_CELL 0.0 end_CELL start_CELL 0.233516 end_CELL start_CELL 11.5365 end_CELL end_ROW end_ARG ] (26)
TABLE I: Simulation parameters
Name Symbol Value
Semi-major axis of target’s orbit a𝑎aitalic_a 7.738e3 km
Chaser mass mcsubscript𝑚𝑐m_{c}italic_m start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT 1500 kg
Control cost penalty γ𝛾\gammaitalic_γ 5
Slack variable penalty ψ𝜓\psiitalic_ψ 750
Thrust bound Umaxsubscript𝑈𝑚𝑎𝑥U_{max}italic_U start_POSTSUBSCRIPT italic_m italic_a italic_x end_POSTSUBSCRIPT 100 N
Position bound rmaxsubscript𝑟𝑚𝑎𝑥r_{max}italic_r start_POSTSUBSCRIPT italic_m italic_a italic_x end_POSTSUBSCRIPT 100 m
Velocity bound vmaxsubscript𝑣𝑚𝑎𝑥v_{max}italic_v start_POSTSUBSCRIPT italic_m italic_a italic_x end_POSTSUBSCRIPT 1.5 m/s
Angular velocity bound ωmaxsubscript𝜔𝑚𝑎𝑥\omega_{max}italic_ω start_POSTSUBSCRIPT italic_m italic_a italic_x end_POSTSUBSCRIPT 0.2 rad/s
Terminal position tolerance ϵpsubscriptitalic-ϵ𝑝\epsilon_{p}italic_ϵ start_POSTSUBSCRIPT italic_p end_POSTSUBSCRIPT 0.35 m
Terminal velocity tolerance ϵvsubscriptitalic-ϵ𝑣\epsilon_{v}italic_ϵ start_POSTSUBSCRIPT italic_v end_POSTSUBSCRIPT 0.03 m/s
Number of docking constraints Nfovsubscript𝑁𝑓𝑜𝑣N_{fov}italic_N start_POSTSUBSCRIPT italic_f italic_o italic_v end_POSTSUBSCRIPT 5
Docking cone half-angle θdocksubscript𝜃𝑑𝑜𝑐𝑘\theta_{dock}italic_θ start_POSTSUBSCRIPT italic_d italic_o italic_c italic_k end_POSTSUBSCRIPT 30
Conservative collision factor αminsubscript𝛼𝑚𝑖𝑛\alpha_{min}italic_α start_POSTSUBSCRIPT italic_m italic_i italic_n end_POSTSUBSCRIPT 1.3
Chaser capture vector in {C}𝐶\{C\}{ italic_C } Dcsubscript𝐷𝑐D_{c}italic_D start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT [002.7]matrix002.7\begin{bmatrix}0&0&2.7\end{bmatrix}[ start_ARG start_ROW start_CELL 0 end_CELL start_CELL 0 end_CELL start_CELL 2.7 end_CELL end_ROW end_ARG ] m
Target capture vector in {T}𝑇\{T\}{ italic_T } Dtsubscript𝐷𝑡D_{t}italic_D start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT [002.7]matrix002.7\begin{bmatrix}0&0&2.7\end{bmatrix}[ start_ARG start_ROW start_CELL 0 end_CELL start_CELL 0 end_CELL start_CELL 2.7 end_CELL end_ROW end_ARG ] m
Safety ellipsoid radial minimum A0min𝐴subscript0𝑚𝑖𝑛A0_{min}italic_A 0 start_POSTSUBSCRIPT italic_m italic_i italic_n end_POSTSUBSCRIPT 15 m
Safety ellipsoid radial maximum A0max𝐴subscript0𝑚𝑎𝑥A0_{max}italic_A 0 start_POSTSUBSCRIPT italic_m italic_a italic_x end_POSTSUBSCRIPT 25 m
Safety ellipsoid cross-track minimum B0min𝐵subscript0𝑚𝑖𝑛B0_{min}italic_B 0 start_POSTSUBSCRIPT italic_m italic_i italic_n end_POSTSUBSCRIPT 10 m
Safety ellipsoid cross-track maximum B0max𝐵subscript0𝑚𝑎𝑥B0_{max}italic_B 0 start_POSTSUBSCRIPT italic_m italic_a italic_x end_POSTSUBSCRIPT 25 m
Maximum iterations imaxsubscript𝑖𝑚𝑎𝑥i_{max}italic_i start_POSTSUBSCRIPT italic_m italic_a italic_x end_POSTSUBSCRIPT 15
Timestep ΔtΔ𝑡\Delta troman_Δ italic_t 1s

VI-A Sequential Convex Procedure Validation

We illustrate the SCP corrective behavior on a single example (Fig. 5). The algorithm slightly adjusts the tail to make it feasible. We observed experimentally the nominal maneuver to require only two burns. Fig. 6 displays the additional burns along the thrust profile required to avoid the large rotating keep-out regions defined by the convex hull of the target. Fig. 7 shows the α𝛼\alphaitalic_α histories; the end of the original trajectory is corrected above the collision threshold after four SCP iterations.

Refer to caption
Figure 6: Nominal thrust trajectory without collision-avoidance constraints follows a distinct two-stage pattern: 1. a de-orbiting burn from the safety ellipse 2. final corrective maneuvers to align the chaser with the target. The sequential convex procedure adds additional ΔVΔ𝑉\Delta Vroman_Δ italic_V maneuvers to the thrust profile to maintain safety.
Refer to caption
Figure 7: The initial convex solution is unsafe towards the end of the maneuver. The SCP corrects the trajectory into the feasible region.

VI-B Convergence and Robustness Analysis

To evaluate the computational and robustness properties of our algorithm, we sampled 250 random attitudes and angular velocity vectors with initial tumble rates up to 10/s{}^{\circ}/\mathrm{s}start_FLOATSUPERSCRIPT ∘ end_FLOATSUPERSCRIPT / roman_s and paired them with 250 random relative orbits with a distance spanning from 15m to 50m. For each sample problem, we incremented the number of timesteps N𝑁Nitalic_N from Nminsubscript𝑁𝑚𝑖𝑛N_{min}italic_N start_POSTSUBSCRIPT italic_m italic_i italic_n end_POSTSUBSCRIPT to Nmaxsubscript𝑁𝑚𝑎𝑥N_{max}italic_N start_POSTSUBSCRIPT italic_m italic_a italic_x end_POSTSUBSCRIPT with a fixed timestep ΔtΔ𝑡\Delta troman_Δ italic_t until the algorithm reached a solution. An infeasible problem was reported if no N𝑁Nitalic_N led to a successful trajectory.

Refer to caption
Figure 8: 250 pairs of initial target tumble rate (up to 10/s{}^{\circ}/sstart_FLOATSUPERSCRIPT ∘ end_FLOATSUPERSCRIPT / italic_s) and attitude and chaser relative orbits (15m to 50m) were sampled as initial conditions. For 93.2% (233) trajectories, the output of the algorithm resulted in a successful and safe capture of the target.

Fig. 8 shows that 93.2% (233) of sampled initial conditions were successful while only 6.8% (17) were infeasible within the bounds of the N𝑁Nitalic_N-search (Nmin=40,Nmax=350formulae-sequencesubscript𝑁𝑚𝑖𝑛40subscript𝑁𝑚𝑎𝑥350N_{min}=40,N_{max}=350italic_N start_POSTSUBSCRIPT italic_m italic_i italic_n end_POSTSUBSCRIPT = 40 , italic_N start_POSTSUBSCRIPT italic_m italic_a italic_x end_POSTSUBSCRIPT = 350). A per-iteration analysis of our algorithm showed an average runtime of 0.15 s per iteration. 8.6% of successes resulted in safe trajectories after only Problem 2, not requiring any corrective actions from the SCP and ensuring global optimality (minimum-fuel). 54.9% of sample trajectories required one SCP iteration for safety, while a cumulative 90.5% of trajectories were safe after five SCP iterations. For a single SCP iteration, the average cumulative runtime was 0.336s, for five SCP iterations it was 0.889s. Much of the fast convergence can be attributed to the near-optimal initial guess after solving Problem 2. Fig. 9 depicts the αksubscript𝛼𝑘\alpha_{k}italic_α start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT evolution of all successful trajectories, confirming their 100% safety guarantees. We observe that using αksubscript𝛼𝑘\alpha_{k}italic_α start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT as a proxy for distance, the trajectory time correlates with the distance which could be used to inform a final time tfsubscript𝑡𝑓t_{f}italic_t start_POSTSUBSCRIPT italic_f end_POSTSUBSCRIPT search online. Note that our method applies in principle to tumble rates ωtmin(ωmax,vmax/(Dc+Dt))delimited-∥∥subscript𝜔𝑡subscript𝜔𝑚𝑎𝑥subscript𝑣𝑚𝑎𝑥subscript𝐷𝑐subscript𝐷𝑡\lVert{}\omega_{t}\rVert{}\leq\min(\omega_{max},v_{max}/(D_{c}+D_{t}))∥ italic_ω start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT ∥ ≤ roman_min ( italic_ω start_POSTSUBSCRIPT italic_m italic_a italic_x end_POSTSUBSCRIPT , italic_v start_POSTSUBSCRIPT italic_m italic_a italic_x end_POSTSUBSCRIPT / ( italic_D start_POSTSUBSCRIPT italic_c end_POSTSUBSCRIPT + italic_D start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT ) ), as faster rates would require the chaser to violate its velocity bounds to achieve soft capture (0.2 rad/s using Table I values).

Refer to caption
Figure 9: Depiction of all α(t)𝛼𝑡\alpha(t)italic_α ( italic_t ) trajectories satisfying the hard collision constraint αk>1subscript𝛼𝑘1\alpha_{k}>1italic_α start_POSTSUBSCRIPT italic_k end_POSTSUBSCRIPT > 1 where αmin=1.3subscript𝛼𝑚𝑖𝑛1.3\alpha_{min}=1.3italic_α start_POSTSUBSCRIPT italic_m italic_i italic_n end_POSTSUBSCRIPT = 1.3 defines a soft safety buffer zone.

VII Discussion and Conclusions

We introduced a fast algorithm to generate safe trajectories for the on-orbit soft capture of tumbling targets with operational and safety constraints. We convexified the field-of-view constraint and employed a sequential convex programming solver to handle the time-varying collision-avoidance constraints. We demonstrated the algorithm’s computational and robustness properties for tumble rates up to 10/s{}^{\circ}{}/\mathrm{s}start_FLOATSUPERSCRIPT ∘ end_FLOATSUPERSCRIPT / roman_s, showcasing its fast convergence and suitability for flight software implementation. Although we decoupled the translational and attitude dynamics, our algorithm still provides feasible reference attitudes for a low-level attitude controller (as a position-dependant function). Furthermore, it quickly detects if the initial conditions are infeasible.

We observed that the infeasible cases tended to occur at both very low and very high angular rates. On the one hand, we were able to render some of these cases feasible by adjusting the discretization parameters (ΔtΔ𝑡\Delta troman_Δ italic_t and N𝑁Nitalic_N). This suggests a need for a free-final time formulation whose naive implementation is, unfortunately, nonconvex. On the other hand, the initial solution obtained by solving Problem 2 dictates the homotopy class of solutions accessible to the SCP part of our algorithm. Given the inherent local nature of the SCP procedure, we suspect that the optimal solution, particularly in low angular velocity regimes, resides in a different homotopy class and remains unreachable with our current method.

Looking ahead, future work includes further investigation of the above limitations, i.e. developing a free-final-time formulation to expand the set of feasible solutions and addressing the homotopy class restriction imposed by our algorithm. On the safety side, plume impingement remains a concern. The plume effects define keep-out zones largely dependent on the thruster configuration used to generate the reference accelerations. Due to its similarity to the collision-avoidance constraints, a similar treatment as in this paper could be envisioned. Finally, uncertainties in the target’s estimated state and dynamics could be accounted for in a variety of ways. An open-source implementation of the algorithm, along with examples, is available at: https://github.com/RoboticExplorationLab/SoftCapture

References

  • [1] W. Fehse, Automated Rendezvous and Docking of Spacecraft.   Cambridge University Press, Nov. 2003.
  • [2] D. C. Woffinden and D. K. Geller, “Navigating the Road to Autonomous Orbital Rendezvous,” Journal of Spacecraft and Rockets, vol. 44, no. 4, pp. 898–909, Jul. 2007, publisher: American Institute of Aeronautics and Astronautics. [Online]. Available: https://arc.aiaa.org/doi/10.2514/1.30734
  • [3] NASA, “On-Orbit Satellite Servicing Study Project Report,” National Aeronautics and Space Administration, Goddard Space Flight Center, Tech. Rep., Oct. 2010.
  • [4] A. Flores-Abad, L. Zhang, Z. Wei, and O. Ma, “Optimal Capture of a Tumbling Object in Orbit Using a Space Manipulator,” Journal of Intelligent & Robotic Systems, vol. 86, no. 2, pp. 199–211, May 2017. [Online]. Available: https://doi.org/10.1007/s10846-016-0417-1
  • [5] J. Virgili-Llop and M. Romano, “Simultaneous Capture and Detumble of a Resident Space Object by a Free-Flying Spacecraft-Manipulator System,” Frontiers in Robotics and AI, vol. 6, 2019.
  • [6] D. C. Sternberg and D. Miller, “Parameterization of Fuel-Optimal Synchronous Approach Trajectories to Tumbling Targets,” Frontiers in Robotics and AI, vol. 5, 2018. [Online]. Available: https://www.frontiersin.org/articles/10.3389/frobt.2018.00033
  • [7] S. Stoneman and R. Lampariello, “A Nonlinear Optimization Method to Provide Real-Time Feasible Reference Trajectories to Approach a Tumbling Target Satellite,” Jun. 2016.
  • [8] K. Albee, C. Oestreich, C. Specht, A. Terán Espinoza, J. Todd, I. Hokaj, R. Lampariello, and R. Linares, “A Robust Observation, Planning, and Control Pipeline for Autonomous Rendezvous with Tumbling Targets,” Frontiers in Robotics and AI, vol. 8, 2021. [Online]. Available: https://www.frontiersin.org/articles/10.3389/frobt.2021.641338
  • [9] D. Malyuta, T. P. Reynolds, M. Szmuk, B. Acikmese, and M. Mesbahi, “Fast Trajectory Optimization via Successive Convexification for Spacecraft Rendezvous with Integer Constraints,” Jun. 2019, arXiv:1906.04857 [math]. [Online]. Available: http://arxiv.org/abs/1906.04857
  • [10] Y. Kim, M. Mesbahi, G. Singh, and F. Hadaegh, “On the Convex Parameterization of Constrained Spacecraft Reorientation,” IEEE Transactions on Aerospace and Electronic Systems, vol. 46, pp. 1097–1109, Aug. 2010.
  • [11] Y.-H. Wu, X.-B. Cao, Y.-J. Xing, P.-F. Zheng, and S.-J. Zhang, “Relative motion coupled control for formation flying spacecraft via convex optimization,” Aerospace Science and Technology, vol. 14, no. 6, pp. 415–428, Sep. 2010. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S1270963810000477
  • [12] X. Shen, S. Diamond, Y. Gu, and S. Boyd, “Disciplined convex-concave programming,” in 2016 IEEE 55th Conference on Decision and Control (CDC), Dec. 2016, pp. 1009–1014. [Online]. Available: https://ieeexplore.ieee.org/document/7798400/similar#similar
  • [13] T. Reynolds and M. Mesbahi, “Small body precision landing via convex model predictive control,” in AIAA SPACE and astronautics forum and exposition, 2017, p. 5179.
  • [14] B. Açıkmeşe, J. Carson, and L. Blackmore, “Lossless Convexification of Nonconvex Control Bound and Pointing Constraints in the Soft Landing Optimal Control Problem,” Control Systems Technology, IEEE Transactions on, vol. PP, pp. 1–1, Nov. 2013.
  • [15] S. Boyd and L. Vandenberghe, Convex optimization.   Cambridge university press, 2004.
  • [16] B. E. Jackson, K. Tracy, and Z. Manchester, “Planning With Attitude,” IEEE Robotics and Automation Letters, vol. 6, no. 3, pp. 5658–5664, Jul. 2021.
  • [17] K. Tracy, T. A. Howell, and Z. Manchester, “Differentiable Collision Detection for a Set of Convex Primitives,” May 2023. [Online]. Available: http://arxiv.org/abs/2207.00669
  • [18] W. H. CLOHESSY and R. S. WILTSHIRE, “Terminal guidance system for satellite rendezvous,” Journal of the Aerospace Sciences, vol. 27, no. 9, pp. 653–658, 1960. [Online]. Available: https://doi.org/10.2514/8.8704
  • [19] NASA, “On-Orbit Satellite Servicing Study Project Report,” National Aeronautics and Space Administration, Goddard Space Flight Center, , Oct. 2010.
  • [20] I. Ross, “How to Find Minimum-Fuel Controllers,” in AIAA Guidance, Navigation, and Control Conference and Exhibit, ser. Guidance, Navigation, and Control and Co-located Conferences.   American Institute of Aeronautics and Astronautics, Aug. 2004. [Online]. Available: https://arc.aiaa.org/doi/10.2514/6.2004-5346
  • [21] S. L. Cleach and Z. Manchester, “Fast solution of optimal control problems with l1 cost,” in Proceedings of AAS/AIAA Astrodynamics Specialist Conference, August 2019.
  • [22] A. Wächter and L. T. Biegler, “On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming,” Mathematical Programming, vol. 106, no. 1, pp. 25–57, Mar. 2006. [Online]. Available: https://doi.org/10.1007/s10107-004-0559-y
  • [23] P. E. Gill, W. Murray, and M. A. Saunders, “SNOPT: An SQP Algorithm for Large-Scale Constrained Optimization,” SIAM Journal on Optimization, vol. 12, no. 4, pp. 979–1006, Jan. 2002. [Online]. Available: https://epubs.siam.org/doi/10.1137/S1052623499350013
  • [24] J. Nocedal and S. J. Wright, Numerical Optimization, 2nd ed.   New York, NY, USA: Springer, 2006.
  • [25] M. ApS, MOSEK Optimizer API for Julia 10.1.24, 2024. [Online]. Available: https://docs.mosek.com/latest/juliaapi/index.html