Delay-Robust Tube Model-Predictive Control for Six-DOF Spacecraft Rendezvous and Docking

Delay-Robust Tube Model-Predictive Control for Six-DOF Spacecraft Rendezvous and Docking Rosalin M. Fernandes Abstract A delay-aware Tube Model Predictive Control (Tube-MPC) architecture is developed, bespoke to the six-degree-of-freedom (6-DOF) dynamics of a chaser spacecraft executing mid-field rendezvous and docking (RVD) with a target vehicle. Tube-MPC is motivated for delay-laden actuator channels, and spacecraft-specific linearization and coupling assumptions are presented. A Gramianbased weight-selection strategy is derived, invariant-set via a Lyapunov ellipsoid is established, and the approach is validated numerically. Results demonstrate that the controlled state trajectory converges to the docking equilibrium while the deviation remains inside an invariant tube despite time-varying actuator delays. Keywords — Tube-MPC, spacecraft rendezvous, delay compensation, robust control, invariant sets I Introduction Rendezvous and docking (RVD) between spacecraft demand precise translational and rotational control in the presence of actuator limits, sensing noise, and communication latency. Delays in the actuator or command channels produce a mismatch between computed control commands and actual applied thrusts, degrading tracking performance and endangering constraint satisfaction. Historically, remedies have included predictorbased compensators, Smith predictors for constant delays, and adaptive or robust MPC variants. However, these methods either use precise delay models or introduce excessive complexity for marginal robustness improvement. Tube-MPC provides an alternative in which an undisturbed nominal trajectory is computed using a standard MPC, and a state-feedback policy is enforced around it to confine the actual trajectory to a bounded tube. This decomposition yields invariance and feasibility guarantees by leveraging a fixed nominal linear model for tube construction while allowing on-line re-linearization for prediction accuracy. Tube-MPC is particularly suitable for spacecraft RVD with bounded, variable actuator delays, offering quantifiable trade-offs between robustness, nominal precision, and computational cost. The key contributions of this study are: •A delay-aware Tube-MPC architecture explicitly tailored to coupled 6-DOF spacecraft dynamics with translational–attitude coupling. •A Gramian-informed weight-normalization framework ensuring balanced modal controllability across heterogeneous state groups. •A Lyapunov-based ellipsoidal invariant certificate validated numerically, demonstrating bounded robustness under variable actuator delays. II Problem statement and methodology A chaser spacecraft operating at mid-field separation from its target is considered, executing an approach and docking manoeuvre within the local-vertical local-horizontal (LVLH) frame. The control objectives are: •to stabilize the relative x-position at 10 m and other positions and attitudes to the docking equilibrium, xk→0which is considered the operating point (x0, u0); •to respect actuator force/torque limits and state constraints; •to guarantee robust invariance to delays modelled as an additive disturbance. The adopted methodology is summarized as follows: 1) The nonlinear 6-DOF plant is linearized about the operating point (x0, u0)to obtain the nominal pair (A, B). 2) An LQR feedback gain K(from the discrete Riccati equation) is computed, and a robust invariant set Eis constructed for the error dynamics. 3) A tightened nominal MPC problem (constraints tightened by E) is solved online, and the combined control u=unom +K(x− xnom)is applied. 4) The model is re-linearized for prediction accuracy, while (Anom, Bnom)is retained for tube verification. III Spacecraft model III-A Nonlinear 6-DOF dynamics The chaser spacecraft is modelled as a rigid body. Using body-frame forces and torques, the continuous nonlinear dynamics are expressed as: ˙p=v, (1) ˙v=1 mR(Θ) Fbody,(2) ˙ Θ=T(Θ) ω, (3) ˙ω=J−1τcmd +τcoupling(Fbody )−ω×Jω, (4) where p, v ∈R3denote position and velocity in LVLH coordinates, Θ = (ϕ, θ, ψ)represents Euler angles, ω∈R3the body rates, R(Θ) the body-to-inertial rotation matrix, T(Θ) the Euler kinematic mapping, mthe spacecraft mass, Jthe inertia matrix, Fbody the commanded body forces, and τcmd the commanded torques. The implemented coupling is realized as: 1) translational acceleration via body-toinertial mapping: acc =1 mR(Θ)Fbody; 2) torque coupling from thruster offset roff: τcoupling =roff ×Fbody; 3) total torque given by τcmd +τcoupling, with angular acceleration J−1(total torque −ω× Jω). This formulation confirms translation–attitude coupling: thrust direction depends on attitude (through R(Θ)), while forces generate moments via thruster offsets. III-B Discrete-time linearization (nominal model) The nonlinear spacecraft model was linearized numerically about the mid-field operating point x0=−50,40.3,−10.5,0.01,0.02,0.04, 0.03,−0.06,0.01,0.005,0.003,0.002⊤ , u0=06. corresponding to a 10mseparation along the xaxis, zero relative velocity, and zero attitude and angular rates. Finite-difference Jacobians of the integrated nonlinear dynamics were evaluated for a sampling period of Ts= 2 s, yielding the discrete model xk+1 =Anomxk+Bnomuk, Anom ∈R12×12, Bnom ∈R12×6. Both matrices are dense due to translational–attitude coupling introduced by the rotation mapping R(Θ) and the thruster-offset moment roff ×Fbody. Their partitioned structure is conveniently expressed as Anom =     App Apv ApΘApω Avp Avv AvΘAvω AΘpAΘvAΘΘ AΘω Aωp Aωv AωΘAωω      , Bnom =     Bp Bv BΘ Bω      , with each block dimensioned 3×3or 3×6as appropriate. Numerically, the following properties were obtained from the full linearization: •∥App∥ ≈ 1.00,∥Apv∥ ≈ 2.00, confirming discrete integrator behaviour in translation; •∥AΘω∥ ≈ 2.00,∥AΘΘ∥ ≈ 1.00, indicating attitude–rate integration; •cross-coupling blocks such as ApΘand Avω have entries of order 10−3–10−2, reflecting weak but non-negligible translational–rotational coupling; •principal gain magnitudes of Bnom satisfy ∥Bv∥≈Ts/m ≃4×10−3for translational inputs and ∥Bω∥≈TsJ−1≃2×10−2for torques; •The pair (Anom, Bnom)is fully controllable (rank(C) = 12). These matrices constitute the fixed nominal linearization used for Tube-MPC synthesis and invariant-set certification, while re-linearization is performed online for prediction accuracy. IV Control architecture: nominal MPC and Tube-MPC IV-A Nominal finite-horizon MPC The nominal MPC solves, at time k, the following receding-horizon problem for Nsteps: min {unom k+i}Jk:= N−1 X i=0 ℓxnom k+i, unom k+i +∥xnom k+N−xref ∥2 P, (5) s.t. xnom k+i+1 =A xnom k+i+B unom k+i,(6) xnom k+i∈Xnom, unom k+i∈Unom.(7) where Q⪰0,R≻0, and Pis the terminal weight (obtained from the discrete Riccati equation or terminal region design). The nominal feasible sets Xnom and Unom are tightened to account for deviation ek. IV-B Tube control law The applied control combines the nominal command and a stabilizing error feedback: uk=unom k+K(xk−xnom k), where Kis the state-feedback matrix computed via the DARE/LQR solution from (Anom, Bnom, Q, R). The closed-loop error dynamics are: Acl := A−BK, ek+1 =Aclek+wk. IV-C Robust positive invariance (RPI) A set Eis robust positively invariant for the error system if: AclE ⊕ W ⊆ E. Under this property, the actual state satisfies: xk∈xnom k⊕ E,∀k, and with tightened constraints: Xnom =X⊖ E, Unom =U⊖KE, recursive feasibility and robust constraint satisfaction follow, as established in [1]. V Weight selection and Gramianbased normalization The selection of the weight matrix Qis based on the discrete controllability Gramian of the nominal linearization. Let (Anom, Bnom)be the nominal pair. The discrete Lyapunov equation: Wc=AnomWcA⊤ nom +BnomB⊤ nom defines modal controllability energy. A statescaling matrix Sis then constructed as: S= diag 1 pdiag(Wc)+ϵ!, where ϵ>0ensures numerical regularization. Normalized weights are obtained as: Q=αQQbias S⊤S Qbias, where Qbias encodes relative state importance and αQis a global tuning multiplier. This Gramianinformed procedure generalizes Bryson’s rule and ensures balanced scaling for tightly coupled modes. VI Invariant set computation Two numerical strategies were employed: VI-A Box fixed-point iteration A component wise iterative bound: s(k+1) =|Acl|s(k)+ ¯w converges to a box RPI only when ∥ |Acl| ∥∞< 1. For the spacecraft system this condition was not met, owing to strong translational–rotational coupling, making the method overly conservative. Table I: Disturbance bounds ¯wand RPI extents ¯sgrouped by physical subsystem. Subsystem ¯w¯sUnits Translation (px, py, pz)0.051–0.052 0.16 m Velocity (vx, vy, vz)0.051–0.052 6.67 m/s Attitude (ϕ, θ, ψ)0.060–0.083 0.09 deg Angular rate (ωx, ωy, ωz)0.060–0.083 0.74 deg/s VI-B Ellipsoidal Lyapunov set An ellipsoidal RPI set is computed from the discrete Lyapunov equation: P=A⊤ cl PAcl +W, W = diag( ¯w2). The invariant ellipsoid: E={e:e⊤P−1e≤1} is positively invariant under Acl for disturbances w∈ W. The principal semi-axes are sellip,i = pλi(P). Robust positive invariance and tightened constraint formulations follow the construction in [2], extending the original Tube-MPC scheme of [1]. VI-C RPI and Robustness Validation Table I summarizes the steady-state RPI extents obtained from the Lyapunov solution, showing all deviations remain within mission-tolerant limits. VII Numerical implementation details •Sampling period: Ts= 2. •Horizon: N= 45. •R and Q constructed via Gramian normalization and bias; tuning multipliers αQ= 20, αR= 1. •Feedback gain: K= (B⊤PB+R)−1B⊤PA, where Pis obtained from the DARE. •Invariant (tube) radius computed from spectral radius: ρ= max i|λi(Acl)|, tube radius ≈∥w∥∞ 1−ρ≈3.0668. The linearization and controller synthesis were implemented using the discrete model derived in Section 3.2 with a sampling period of Ts= 2 s and horizon N= 45. The main numerical properties of the resulting Tube-MPC setup are summarized below. •State and input dimensions: nx= 12,nu= 6; total decision variables in the quadratic program: 822. •Control authority: umin = [−20,−20,−20,−5,−5,−5]⊤, umax = [20,20,20,5,5,5]⊤. The maximum forces and torques recorded during simulation were |uF|max = 7.58 N and |uτ|max = 0.13 N m, both within actuator limits. •Algebraic Riccati solution: The discretetime Riccati equation yielded Pwith eigenvalue range and condition number as shown in Table II. •Closed-loop stability: The closed-loop matrix Acl =A−BK had eigenvalues grouped as: λ(Acl) = 0.9639 ±0.0182i(×6) 0.8507 ±0.0682i(×6) , ρ(Acl)=0.964 <1. all strictly within the unit circle, confirming asymptotic discrete-time stability. •System conditioning: cond(A) = 5.83 with row norms in (0.999,1.618,2.236), indicating well-scaled dynamics. These diagnostics confirm numerical stability of the Riccati and MPC solvers, bounded control effort, and well-conditioned weighting matrices. VIII Results VIII-A Qualitative description Figure 1a presents translational and rotational state trajectories for variable delay dk∈[0,3]. The nominal trajectory is shown for reference. The actual closed-loop trajectories converge to the docking equilibrium and remain confined within the certified ellipsoidal tube. Table II: Closed-loop and numerical conditioning metrics for Tube-MPC synthesis. Metric Value Interpretation cond(Q) 3.20 ×103Well-scaled state weighting cond(R) 10.0Balanced input penalization λmin(P) 0.0398 Positive definite Lyapunov matrix λmax(P) 528.68 Positive definite Lyapunov matrix ρ(Acl) 0.964 Discrete-time stability margin Decision vars. 822 Quadratic program dimension (a) State trajectories under Tube-MPC: positions (solid) and velocities (dashed) for Translation (left). angular positions (solid) and velocities (dashed) for Rotation(right); states converge to the docking equilibrium. (b) Evolution of deviation norms and invariant bounds. The translational error norm (blue) enters and remains within the bound (red dashed line at 3.0668) after t≈200 s. Rotational deviation (orange) remains within bounds for the entire horizon. VIII-B Invariant containment and entry time Figure 1b illustrates the norms of translational and rotational deviations along with invariant bounds. The translational error norm falls below the invariant radius (approximately 3.07) at t≈200 s and remains within thereafter, while rotational deviations remain an order of magnitude smaller. VIII-C Benchmarking against delay-free nominal MPC Performance can be quantified via RMS error, steady-state error, and control effort for: •nominal delay-free MPC baseline (no tube, no delay); •nominal MPC with delay (no tube); •Tube-MPC under the same delay realization. Figure 2: Baseline MPC without delay used as a benchmark Figure 3: Nominal MPC under a two-step actuator delay without tube constraints. The response remains stable but exhibits significant angular overshoot (≈0.10×rad) due to delay-induced phase lag. The percentage improvement is 11.1% and computed as: ImprovementRMS = 100%·RMSbaseline −RMStube RMSbaseline . In Fig. 3; however, its transient behaviour exhibits pronounced angular overshoots, with peak roll and pitch excursions exceeding 0.10 rad. These peaks arise from the phase lag between the predicted and realized control actions, a characteristic vulnerability of delay-afflicted feedback systems. In contrast, the Tube-MPC maintains comparable steady-state precision while reducing these transient peaks by 10–15% and completely suppressing oscillatory rebound in the angular channels. This demonstrates that the Tube formulation not only preserves stability under delay but also mitigates delay-induced overshoot through its invariant feedback structure. IX Discussion: trade-offs and limitations •Robustness vs nominal performance: TubeMPC sacrifices nominal optimality (tighter constraint sets, feedback correction) for guaranteed robustness under delay; this may be quantified through RMS and convergencetime metrics. •Fixed nominal linearization: Certification uses (Anom, Bnom). Re-linearization within the MPC improves predictions but invalidates invariance proofs if used to rebuild the tube online. This represents a deliberate and defensible design trade-off. •Computational cost: Tube-MPC requires online solution of the nominal MPC, with complexity comparable to standard MPC using tightened constraints. Actual solve-time statistics are reported in the results. Table III: Improvement in attitude–tracking performance under Tube-MPC (delay-aware) relative to nominal MPC. Errors(rad) Nominal Tube-MPC Impro.(%) RMS roll 0.0306 0.0234 23.3 RMS pitch 0.0258 0.0198 23.4 RMS yaw 0.0189 0.0171 9.2 Max roll 0.0761 0.0647 15.0 Max pitch 0.0667 0.0600 10.0 Max yaw 0.0521 0.0511 1.9 Control-effort trade-off: Metric Nominal Tube-MPC RMS control effort 1.4157 2.2030 Peak control effort 6.3394 10.5778 Tube-MPC consumes greater control energy and peak actuation (≈1.5×nominal) due to the additional feedback action that maintains invariance under delay perturbations. The Tube-MPC achieved consistent improvements in attitude-tracking accuracy (Table III), reducing RMS roll and pitch errors by approximately 23% and yaw RMS by about 9%, while also lowering peak deviations across all axes. Steady-state offsets were effectively eliminated for both controllers. The enhanced performance is obtained at the expense of higher control usage, as summarized in the trade-off table above, where RMS and peak control efforts increase by roughly 50–65%. This behaviour exemplifies the characteristic robustness–energy exchange of Tube-MPC: additional actuation is expended to guarantee constraint satisfaction and delay-tolerant invariance. X Conclusions The formulation is directly extensible to delayafflicted formation control, proximity operations, and high-precision pointing systems, where actuator latency is an unavoidable performance constraint. Simulation results confirm that the actual trajectory approaches the nominal trajectory and remains inside the certified invariant tube despite bounded, time-varying actuator delays. Future work includes refinement of the ellipsoidal certificate using less conservative disturbance models and validation using Monte Carlo methods. The developed formulation provides a computationally tractable route for embedding robust delay compensation in on-board MPC architectures for spacecraft and robotic systems. References [1] D. Q. Mayne, M. M. Seron, and S. V. Rakovic, “Robust model predictive control of constrained linear systems with bounded disturbances,” Automatica, vol. 41, no. 2, pp. 219–224, 2005. [2] S. V. Rakovic, W. S. Levine, and D. Q. Mayne, “Model predictive control with robustness and constraint satisfaction,” Springer Handbook of Model Predictive Control, pp. 639–662, 2012.