SmoothedImpulsiveBurn#
- class SmoothedImpulsiveBurn(spacecraft: Body, maneuvers, base_frame: Frame = None, t0: float = None, smoothing_factor: float = 1.0)#
Bases:
DynamicModelModels impulsive maneuvers as continuous tanh-shaped acceleration pulses.
A true impulsive maneuver applies an instantaneous velocity increment at \(t_{\Delta v}\), which makes the state discontinuous and therefore forces the trajectory to be split into legs around every burn (see
MissionSequence). This model removes that constraint by replacing the Dirac impulse with a narrow, continuously differentiable pulse, so a single uninterrupted propagation can carry the maneuver and the variational equations can be integrated straight through it.The pulse is built from the product of two opposite-signed hyperbolic tangents, which yields a localized peak that decays to zero away from the maneuver epoch:
\[s_{\Delta v}(t) = -\tanh^2\!\big(a (t - t_{\Delta v})\big) + 1 = \operatorname{sech}^2\!\big(a (t - t_{\Delta v})\big)\]where \(a\) is the smoothing factor. Integrating from the start of the arc gives the normalization
\[\int_{t_0}^{\infty} s_{\Delta v}\,\mathrm{d}t = \frac{1}{a}\Big(1 + \tanh\!\big(a (t_{\Delta v} - t_0)\big)\Big)\]so that dividing by it makes the pulse deliver exactly \(\Delta v\):
\[\mathbf{u}^{\text{smooth}}_{\Delta v}(t) = \frac{-\tanh^2\!\big(a (t - t_{\Delta v})\big) + 1} {\frac{1}{a}\big(1 + \tanh(a (t_{\Delta v} - t_0))\big)}\, \Delta \mathbf{v}\]As \(a \to \infty\) this approaches \(\Delta\mathbf{v}\,\delta(t - t_{\Delta v})\) in the distributional sense. In practice \(a\) must stay finite: too large a value produces gradients steep enough to collapse the integrator step size. See
smoothing_factorfor guidance on choosing it.The extended dynamics augment the orbital state with the maneuver parameters \(\boldsymbol{\theta}_{\Delta v} = [\Delta\mathbf{v},\, \delta t_{\Delta v}]\), whose own dynamics are zero:
\[\begin{split}\dot{\mathbf{x}}_{\Delta v} = \begin{bmatrix} \dot{\mathbf{X}} \\ \dot{\boldsymbol{\theta}}_{\Delta v} \end{bmatrix} = \begin{bmatrix} \mathbf{f} \\ \mathbf{0}_{4\times 1} \end{bmatrix} + \mathbf{B}_{\Delta v}\, \mathbf{u}^{\text{smooth}}_{\Delta v}\end{split}\]with \(\mathbf{B}_{\Delta v} = [\mathbf{0}_{3\times3};\, \mathbf{I}_3;\, \mathbf{0}_{4\times3}]\). Because this model returns an acceleration, the propagator supplies \(\mathbf{B}_{\Delta v}\) implicitly and the maneuver parameters are static state entries.
- Parameters:
- spacecraftBody
The maneuvering body. Its
spice_idkeys the parameter lookup.- maneuversImpulsiveBurn or list of ImpulsiveBurn
The maneuvers to smooth. Each supplies the nominal delta-v vector and nominal epoch
t_burn(TDB seconds past J2000). These are the same objects a leg-basedMissionSequencewould consume.- base_frameFrame, optional
Reference frame the delta-v vectors are expressed in. Defaults to J2000.
- t0float, optional
Arc start epoch (TDB seconds past J2000) used as the lower limit of the normalization integral. When
None, the pulse is normalized over the whole real line (equivalent to assuming the maneuver sits many pulse widths after the arc start), which is the usual case.- smoothing_factorfloat, optional
The smoothing factor \(a\), in \(s^{-1}\). Defaults to
1.0.
- Attributes:
base_frameReference frame the delta-v vectors are expressed in.
estimable_parameter_namesThe solve-for parameter names this model responds to.
maneuversThe list of
ImpulsiveBurnobjects being smoothed.smoothing_factorThe smoothing factor \(a\), in \(s^{-1}\).
spacecraftThe maneuvering body.
t0Arc start epoch used as the lower limit of the normalization integral.
Methods
compute_acceleration(position, ...)Computes the smoothed acceleration from all modeled maneuvers.
criterion_residual(a, dt_dv[, eta, t_dv, t0])The criterion function \(g_{\Delta v}\), Eq.
mapped_dv_ratio(a, dt_dv[, t_dv, t0])Mapped velocity deviation at the interval midpoint, per delta-v, Eq.
optimal_smoothing_factor(dt_dv[, eta, t_dv, ...])Solve \(g_{\Delta v} = 0\) for the smoothing factor.
total_delivered_dv(epochs[, current_state, body])Numerically integrate the pulse over
epochsto check delivered delta-v.- Raises:
- SmoothedImpulsiveBurn.NoManeuversError
If
maneuversis empty.- SmoothedImpulsiveBurn.InvalidSmoothingFactorError
If
smoothing_factoris not strictly positive.
See also
scarabaeus.ImpulsiveBurnInstantaneous delta-v applied across legs.
scarabaeus.FiniteBurnContinuous finite-duration thrust model.
Notes
For maneuver
i(1-based), this model reads two solve-for parameters from the propagated state, falling back to the values supplied at construction when they are absent:dv_smooth{i}The delta-v vector (size 3, units km/s).
dt_dv{i}A correction to the nominal maneuver epoch (size 1, units sec), so the effective epoch is
t_burn + dt_dv. Estimating an offset rather than the absolute epoch keeps the filter well conditioned: an absolute TDB epoch is order \(10^{8}\) s, and its correction is typically order seconds. The partial is unchanged by this shift, since \(\partial t_{\Delta v} / \partial \delta t_{\Delta v} = 1\).
When exactly one maneuver is modeled, the unsuffixed names
dv_smoothanddt_dvare also accepted.References
[1]Tanh-based smoothing parameterization of impulsive maneuvers for state-transition-tensor propagation, Eqs. (17)-(20).
- classmethod criterion_residual(a: float, dt_dv: float, eta: float = 0.9, t_dv: float = None, t0: float = None) float#
The criterion function \(g_{\Delta v}\), Eq. (29).
Requires the mapped velocity deviation at the midpoint to be a fixed fraction of the maneuver, with the sign that matches the real deviation: \(\delta v_x(t') = -\eta \Delta v_x\) for \(\delta t_{\Delta v} > 0\) and \(+\eta \Delta v_x\) for \(\delta t_{\Delta v} < 0\). The residual is
mapped_dv_ratio()plus \(\eta\,\operatorname{sign}(\delta t_{\Delta v})\).- Parameters:
- afloat
Smoothing factor, in \(s^{-1}\).
- dt_dvfloat
Maneuver-epoch deviation \(\delta t_{\Delta v}\), in seconds.
- etafloat, optional
Target fraction \(\eta\) of the maneuver. Defaults to
0.9.- t_dv, t0float, optional
Maneuver epoch and arc start, passed to
mapped_dv_ratio(). Defaults toNone.
- Returns:
- float
Zero when
asatisfies the midpoint relation.
- classmethod optimal_smoothing_factor(dt_dv: float, eta: float = 0.9, t_dv: float = None, t0: float = None, a_bounds: tuple = (0.001, 4.0), n_grid: int = 4000, branch: str = 'low') float#
Solve \(g_{\Delta v} = 0\) for the smoothing factor.
Implements the grid search for root approximation (GSRA):
ais discretized overa_bounds,criterion_residual()is evaluated on the grid, and each sign change is refined by bisection. No analytic solution exists.In the \(S_2 \to 1\) limit the criterion collapses to a function of the single group \(u = a\,\delta t_{\Delta v}/2\),
\[u \operatorname{sech}^2(u)\,(1 + 2u\tanh u) = \eta\]so the smoothing factor scales as \(a \propto 1/\delta t_{\Delta v}\) and is independent of the maneuver magnitude, the maneuver epoch, and the orbit. Because \(h(u)\) rises then falls, the equation has two roots, selected by
branch.- Parameters:
- dt_dvfloat
The maneuver-epoch uncertainty to tune for, in seconds. Sign is irrelevant: the criterion is symmetric, so
|dt_dv|is used.- etafloat, optional
Adjustment coefficient. Defaults to
0.9, reported as the value minimizing the position-deviation residual. Must be belowETA_MAX.- t_dv, t0float, optional
Passed through to
mapped_dv_ratio(); see the note there on why they barely matter.- a_boundstuple of float, optional
(a_min, a_max)search interval, defaults to(1e-3, 4.0).- n_gridint, optional
Number of grid points for the search. Defaults to
4000.- branch{‘low’, ‘high’}, optional
Which of the two roots to return.
'low'(default) gives the wider pulse, which is the numerically safer choice: it stays further from the step-resolution limit described insmoothing_factor.
- Returns:
- float
The smoothing factor, in \(s^{-1}\).
- Raises:
- SmoothedImpulsiveBurn.NoSmoothingFactorRootError
If
etais unreachable or no root lies insidea_bounds.
Examples
>>> SmoothedImpulsiveBurn.optimal_smoothing_factor(100.0) 0.015526...
- static mapped_dv_ratio(a: float, dt_dv: float, t_dv: float = None, t0: float = None) float#
Mapped velocity deviation at the interval midpoint, per delta-v, Eq. (28).
When only a maneuver-epoch deviation \(\delta t_{\Delta v}\) is present, the mapped velocity deviation follows from the extended first- and second-order sensitivities
\[\delta \dot v_x \approx A^{\Delta v}_{4,10}\,\delta t_{\Delta v} + \tfrac{1}{2}\Lambda^{\Delta v}_{4,10,10}\,\delta t_{\Delta v}^2\]with \(A^{\Delta v}_{4,10} = \Delta v_x\,\partial_{t_{\Delta v}} P\) and \(\Lambda^{\Delta v}_{4,10,10} = \Delta v_x\, \partial^2_{t_{\Delta v}} P\), where \(P\) is this model’s normalized pulse integral. Carrying those derivatives through and integrating in time gives Eq. (28),
\[\frac{\delta v_x}{\Delta v_x} = \frac{a}{S_2 + 1} \Big[S_1^2 + (S_2 - 1) S_1 - S_2\Big] \Big(\delta t_{\Delta v} + a S_1 \delta t_{\Delta v}^2\Big)\]with \(S_1 = \tanh(a(t - t_{\Delta v}))\) and \(S_2 = \tanh(a(t_{\Delta v} - t_0))\), evaluated here at the midpoint \(t' = t_{\Delta v} + \tfrac{1}{2}\delta t_{\Delta v}\) so that \(S_1 = \tanh(a\,\delta t_{\Delta v}/2)\).
\(\Delta v_x\) divides out, which is why the criterion is independent of the maneuver magnitude.
- Parameters:
- afloat
Smoothing factor, \(s^{-1}\).
- dt_dvfloat
Maneuver-epoch deviation \(\delta t_{\Delta v}\), in seconds.
- t_dv, t0float, optional
Maneuver epoch and arc start. They enter only through \(S_2\); when either is omitted, \(S_2 = 1\) is used, the limit reached whenever the maneuver sits more than a few pulse widths after the arc start. This is why \(t_{\Delta v}\) has almost no influence on the resulting smoothing factor.
- Returns:
- float
\(\delta v_x(t') / \Delta v_x\), dimensionless.
- compute_acceleration(position: ndarray, current_state: dict, epoch: float, body: Body) ndarray#
Computes the smoothed acceleration from all modeled maneuvers.
- Parameters:
- positionnp.ndarray
3x1 spacecraft position vector. Unused: the pulse depends only on time and the maneuver parameters.
- current_statedict
Dictionary containing the current state parameters.
- epochfloat
Current simulation epoch (TDB seconds past J2000).
- bodyBody
The body whose maneuver parameters are read.
- Returns:
- accnp.ndarray
(3,)acceleration in km/s^2.
- total_delivered_dv(epochs: ndarray, current_state: dict = None, body: Body = None) ndarray#
Numerically integrate the pulse over
epochsto check delivered delta-v.Diagnostic helper for the analytic normalization: on a grid fine enough to resolve the pulse (spacing well below \(1/a\)) this returns the commanded delta-v to roundoff, confirming Eq. (18).
Warning
This does not predict whether a propagation will resolve the pulse. An adaptive integrator chooses its own internal steps, which are unrelated to the epochs a run reports on: a pulse this helper cannot see on a coarse output grid may still be integrated perfectly. The only reliable check that a given
smoothing_factoris safe is to compare a propagation against the equivalent two-leg impulsive answer once, astests/python_tests/unit_testing/dynamics_tests/test_SmoothedImpulsiveBurn.pydoes.- Parameters:
- epochsnp.ndarray
Monotonic grid of epochs (TDB seconds past J2000).
- current_statedict, optional
Parameter overrides; defaults to the nominal values.
- bodyBody, optional
Owner for parameter lookup. Defaults to the model’s spacecraft.
- Returns:
- np.ndarray
(3,)integrated delta-v in km/s.
- ETA_MAX = 1.098686#
Largest attainable
eta. The criterion function \(h(u) = u\\,\\operatorname{sech}^2 u\\,(1 + 2u\\tanh u)\) peaks ath = 1.098686foru = 1.182514, so no smoothing factor can satisfy the midpoint relation beyond this. The published choice, 0.9, sits comfortably inside.
- property estimable_parameter_names: tuple[str, ...]#
The solve-for parameter names this model responds to.
- Returns:
- tuple of str
('dv_smooth1', 'dt_dv1', 'dv_smooth2', ...), plus the unsuffixed aliases when a single maneuver is modeled.
- property maneuvers: list#
The list of
ImpulsiveBurnobjects being smoothed.
- property smoothing_factor: float#
The smoothing factor \(a\), in \(s^{-1}\).
The pulse has an effective width of order \(1/a\) and a peak acceleration of order \(a\,\|\Delta v\|/2\).
Accuracy against the true impulse improves as \(O(1/a^2)\) (the pulse is symmetric, so its leading error term is quadratic), but only while the integrator still resolves it. Past that point the adaptive step controller steps clean over the pulse, its error estimate never spikes, so it never refines, and the maneuver is silently under- or un-applied. That failure is not gradual: in the LEO case in
tests/python_tests/unit_testing/dynamics_tests/test_SmoothedImpulsiveBurn.py, the final position error falls from 1.1e-2 km to 7.2e-4 km to 2.9e-5 km for \(a = 0.05, 0.2, 1.0\,s^{-1}\), then jumps to 1.7e+1 km at \(a = 5\,s^{-1}\).Scarabaeus exposes no maximum-step control, so \(a\) is the only lever: tightening
atol/rtoldoes not rescue an unresolved pulse. Choose the smallest \(a\) that meets the accuracy target, and validate it once against the equivalent two-leg impulsive propagation before trusting a new configuration.Important
Trajectory fidelity and deviation-mapping fidelity pull \(a\) in opposite directions, and
optimal_smoothing_factor()optimizes the second, not the first. It returns \(a \propto 1/\delta t_{\Delta v}\), which for realistic epoch uncertainties is a far wider pulse than pure trajectory reproduction wants. Measured on the LEO case:\(\delta t\)
\(a^*\)
trajectory error
50 s
0.0311
2.97e-02 km
100 s
0.0155
1.18e-01 km
300 s
0.0052
1.04e+00 km
against 2.9e-05 km at \(a = 1\). Use the criterion when the deliverable is a sensitivity mapping over an uncertain maneuver epoch; use a large \(a\) when the deliverable is the trajectory itself. If you need both, they are different runs.