Residual-Aware Control: Bounds and Decisions for Nonlinear Motion
From Taylor Remainders to Task-Relevant Error Budgets
Executive Summary
A local linear model tells us how a small change in state or control changes a motion. Its remainder tells us how much that prediction misses for a specified model, comparison trajectory, coordinate system, and interval. This makes residual analysis useful for deciding where to refine a trajectory or restrict an update. It does not turn every measured error into a sensor of intrinsic geometric curvature.
The important connection is between three quantities: the nonlinear remainder generated now, its amplification or attenuation during the remaining motion, and its effect on the eventual task. In a golf swing, an error that is small in a joint-coordinate norm can still matter greatly for clubface orientation or strike location. A larger deviation in another direction may scarcely affect the shot. In walking, a small error in touchdown timing can change the contact sequence. The norm alone cannot decide which discrepancy deserves attention.
This article develops a propagated flow-error bound, a dimensionally consistent interval rule, and the additional conditions needed for observation, switching, and robust tubes. It includes reproducible scalar calculations and distinguishes them from proposed quadrotor, humanoid, and golf evaluations. It makes no claim that a residual-adaptive solver has demonstrated a particular speedup, fall reduction, or improvement in golfers.
Readers should know state-space models, multivariable Taylor expansion, and basic linear feedback. The series introduction provides the broader context. Work through the two appendices if you want to check the bounds rather than use them conditionally. State-transition matrices describe first variations; Hessians describe how those first variations change. These are different roles.
1 Part I: Theory
1.1 Residual Definition and Quantitative Bounds
1.1.1 Define the Comparison Before Interpreting the Error
Start within one smooth mode of a deterministic model. Use scaled coordinates \(z=S_x^{-1}x\) and \(v=S_u^{-1}u\), with fixed, nonsingular diagonal scale matrices. For example, angles, angular velocities, activation states, and torques need distinct reference scales. Write
\[ \dot z=F(t,z,v),\qquad \dot{\bar z}=F(t,\bar z,\bar v). \]
Time retains its declared physical unit. Euclidean norms of \(z\) and \(v\) are dimensionless; \(F\), its scaled-coordinate Hessian, and its Jacobian have units of inverse time. Choosing scales is a modeling decision, not cosmetic formatting: changing them changes which errors the norm emphasizes.
Let \(e=z-\bar z\) be the actual state deviation and \(\delta v=v-\bar v\) the actual input difference. Define the first-order comparison by
\[ \begin{aligned} \dot\eta &= A(t)\eta+B(t)\delta v(t),\quad \eta(t_0)=e(t_0),\\ A(t)&=D_zF(t,\bar z,\bar v),\quad B(t)=D_vF(t,\bar z,\bar v). \end{aligned} \]
The input difference is prescribed for this comparison. In an offline calculation it can be the same recorded input applied to two model trajectories. If it depends on the actual state, computing this comparison online requires that input history; it is not automatically an independent predictor.
\[ z(t)=\bar z(t)+\eta(t)+r(t),\qquad r=e-\eta. \]
This \(r\) is the flow remainder. It differs from tracking error \(e\), a one-step integration defect, and a measured prediction innovation. For a fixed finite horizon, a uniformly smooth flow, a common domain of existence, and an initial/input perturbation family of size \(\alpha\), the first variation has
\[ \|r(t)\|=O(\alpha^2). \tag{1}\]
The constant depends on the horizon and dynamics. This asymptotic statement neither specifies an acceptable error nor stays uniform as the horizon approaches a singularity.
1.1.2 Hessian-Based Residual Bounds
Assume \(F\) is twice continuously differentiable in \((z,v)\). At each time, the segment joining \((\bar z,\bar v)\) and \((z,v)\) must lie inside the stated domain. Put \(w=(e,\delta v)\). Taylor’s integral formula gives an exact remainder without assuming a third derivative:
\[ \begin{aligned} N(t,w)&=F(t,z,v)-F(t,\bar z,\bar v)-Ae-B\delta v,\\ N(t,w)&=\int_0^1(1-s)D^2_{(z,v)}F(t,\bar z+se,\bar v+s\delta v)[w,w]\,ds. \end{aligned} \tag{2}\]
Here a vector-valued Hessian is a bilinear map. Suppose its induced norm satisfies \(\|D^2F[a,b]\|_2\le M(t)\|a\|_2\|b\|_2\) throughout the joined segment, or throughout a region known to contain it. Then
\[ \|N(t,w)\|_2\le\frac{M(t)}2\bigl(\|e(t)\|_2^2+\|\delta v(t)\|_2^2\bigr). \tag{3}\]
Subtract the comparison dynamics from the actual deviation dynamics:
\[ \dot r=A(t)r+N(t,w),\qquad r(t_0)=0. \tag{4}\]
Let \(\Phi(t,s)\) solve \(\partial_t\Phi=A(t)\Phi\), \(\Phi(s,s)=I\). Variation of constants gives
\[ r(t)=\int_{t_0}^t\Phi(t,s)N(s,w(s))\,ds. \tag{5}\]
Under the smoothness, domain, norm, and comparison assumptions above,
\[ \|r(t)\|_2\le\frac12\int_{t_0}^t \|\Phi(t,s)\|_2M(s) \bigl(\|e(s)\|_2^2+\|\delta v(s)\|_2^2\bigr)\,ds. \tag{6}\]
The deviation inside this bound is the actual deviation. Substituting the linear prediction requires a separate bound on its difference from the actual deviation. Dropping \(\Phi\) requires an established nonexpansive propagator in the chosen norm.
The theorem describes how the model’s local nonlinearity affects a finite motion. It does not cover omitted forces or numerical integration error unless their contributions are introduced explicitly. A measured disturbance term \(d\) would add \(\int\Phi d\) to the error equation, with its own assumptions and bound.
Feedback requires particular care. If \(\delta v=Ke\) with a known constant or time-dependent linear gain, then \(\|w\|^2\le(1+\|K\|^2)\|e\|^2\) for the open-loop comparison just defined. Its forcing bound becomes
\[ \|N\|\le\frac{M}2(1+\|K\|^2)\|e\|^2. \tag{7}\]
Alternatively, define the closed-loop field \(F_K(t,z)=F(t,z,\bar v+K(z-\bar z))\) and linearize that field. Its propagator uses \(A+BK\), and its Hessian must include the feedback substitution. A nonlinear feedback law contributes its own second derivatives. One cannot use a closed-loop propagator with an open-loop remainder definition without rederiving the comparison. The appendix makes this distinction explicit.
1.1.3 Computing the Hessian Bound
For a control-affine field \(F(z,v)=f(z)+\sum_jg_j(z)v_j\),
\[ \begin{aligned} D^2_{zz}F&=D^2f+\sum_jv_jD^2g_j,\\ D^2_{zv_j}F&=Dg_j,\qquad D^2_{vv}F=0. \end{aligned} \tag{8}\]
Thus control-affine does not mean jointly linear in state and input. State-dependent actuation produces mixed derivatives, even though input-input derivatives vanish. With \(F(z,v)=zv\), for instance, the mixed remainder is \(e\delta v\).
If output component \(F_k\) has Hessian matrix \(H_k\), the estimate \(\sqrt{\sum_k\|H_k\|_2^2}\) bounds the bilinear norm for Euclidean inputs and outputs. The largest component norm alone supplies an infinity-output bound; multiplying it by \(\sqrt n\) gives a conservative Euclidean-output bound for \(n\) components. Mixing these conventions silently changes the theorem.
Finite differences and automatic differentiation evaluate derivatives at points. A sample maximum can miss a larger derivative between samples. Certification over a region requires an analytical bound, verified interval calculation, or another justified enclosure. A dense sampled maximum is useful as an empirical indicator, provided the article or controller labels it that way. Differentiate the equations actually used, including coordinate changes, input maps, and smooth actuator dynamics.
Pendulum example. With point mass \(m\), length \(L\), viscous rotational damping coefficient \(b\) in torque per angular velocity, and applied torque \(u\),
\[ \dot\theta=\omega,\qquad \dot\omega=-\frac gL\sin\theta-\frac b{mL^2}\omega+\frac u{mL^2}. \]
The acceleration component’s state Hessian is
\[ H_{\dot\omega}= \begin{bmatrix}(g/L)\sin\theta&0\\0&0\end{bmatrix}. \tag{9}\]
Its sign follows by differentiating \(-\sin\theta\) twice. In unscaled coordinates the acceleration remainder from an angular displacement \(\Delta\theta\) obeys
\[ |N_{\dot\omega}|\le\frac{g}{2L}|\Delta\theta|^2. \tag{10}\]
At \(L=1\,\mathrm m\) and \(|\Delta\theta|\le0.1\) rad, using \(g=9.81\,\mathrm{m\,s^{-2}}\) for this rounded illustration, the acceleration-remainder bound is \(0.04905\,\mathrm{rad\,s^{-2}}\). Integrating a constant acceleration remainder for \(0.1\) s in a double-integrator comparison would give
\[ |r_\omega|\le0.004905\,\mathrm{rad\,s^{-1}},\qquad |r_\theta|\le0.00024525\,\mathrm{rad}. \tag{11}\]
These are different components with different units. For the pendulum’s actual variational dynamics, use its coupled transition matrix instead of the double-integrator calculation. An angle/rate vector cannot be assigned a single unit of radians without scaling.
1.2 A Counterexample That Exposes Missing Propagation
Consider the dimensionless teaching model \(\dot x=x+x^2\), nominal trajectory \(\bar x=0\), and initial deviation \(\epsilon>0\). Before its finite-time pole,
\[ x(t)=\frac{\epsilon e^t}{1-\epsilon(e^t-1)},\qquad \eta(t)=\epsilon e^t. \]
Its Hessian is exactly \(2\), and its transition matrix is \(e^{t-s}\). The tempting expression obtained by dropping propagation and substituting \(\eta\) is \(B_{\mathrm{old}}=\epsilon^2(e^{2T}-1)/2\). But the actual remainder is
\[ r(T)=\frac{\epsilon^2e^T(e^T-1)}{1-\epsilon(e^T-1)}. \]
As \(\epsilon\to0\), \(r/B_{\mathrm{old}}\to2e^T/(e^T+1)>1\) for every \(T>0\). The proposed upper bound fails even for arbitrarily small perturbations. At \(\epsilon=0.01\), \(T=1\), the actual remainder is approximately \(0.000475\), while the old expression is approximately \(0.000319\).
For a simple valid enclosure, the positive solution increases, so \(x(s)\le x(T)\) throughout this interval. The corrected conditional bound is \(x(T)^2(e^T-1)\), approximately \(0.001314\). This is deliberately conservative. Using the actual time-dependent deviation in the propagated integral recovers the exact remainder here because all terms are scalar and positive. The implementation section reproduces the comparison; its numerical tests also compare the closed form with an independently integrated ODE.
2 Residuals as Control Signals
2.1 Real-Time Residual Monitoring
Three comparisons answer three different questions:
| Quantity | Comparison | What It Can Establish |
|---|---|---|
| Model flow remainder | Nonlinear model versus its first variation, same initial/input experiment | Validity of the chosen local approximation |
| Tracking discrepancy | Estimated state versus the desired or nominal state | Departure from the selected motion |
| Prediction innovation | New measurement versus a prediction made before that measurement | Consistency of prediction and observation under an estimation model |
A one-step first-order prediction can be constructed from the previous estimate:
\[ \hat z_{k+1|k}=\bar z_{k+1} +\Phi_k(\hat z_{k|k}-\bar z_k)+\Gamma_k\delta v_k. \tag{12}\]
For piecewise-constant input deviations, \(\Gamma_k\) is the appropriate integrated input sensitivity. If inputs vary within the interval, use their convolution with the transition matrix. The prediction must include the affine offset through the evolving nominal trajectory. It must not reuse the current measurement as its own initial condition.
For measurement function \(h\), the innovation is
\[ \nu_{k+1}=y_{k+1}-h(\hat z_{k+1|k}). \tag{13}\]
It combines state-estimation error, model discrepancy, unmodeled inputs, numerical error, measurement error, and linearization remainder. An exact measurement of every state still does not identify which omitted physical force caused a mismatch. A filtered innovation is useful evidence about consistency; filtering alone cannot make indistinguishable causes identifiable.
For a locally linear observation and independent prior state error and measurement noise, the usual innovation covariance approximation is \(S=HP^-H^T+R\). If those errors are correlated, the corresponding cross-covariance terms must be included with the signs implied by the error convention. A normalized statistic \(\nu^TS^{-1}\nu\) can help compare measurements of different units. A chi-squared interpretation additionally requires the stated distribution and calibration; it is not supplied by the algebra alone.
2.2 Threshold-Based Decision Making
Set thresholds from a task budget and an estimation model. A large discrepancy may call for a better state estimate, a smaller proposed control update, a different contact mode, or a fallback policy. It does not necessarily call for a shorter integration step. If the cause is a missing force, integrating the same wrong model more accurately can leave the prediction just as wrong.
Record which residual triggered an action. A controller report should contain the prediction timestamp, interval, reference, scales, model version, norm, uncertainty estimate, and resulting action. This allows a reviewer to determine whether a threshold crossing preceded the problem or merely described it afterward.
2.3 Example 1.2: Quadrotor Flip
A flip illustrates why angle, rate, and approximation quality must be separated. A large attitude excursion makes a single small-angle approximation unsuitable, but a moving local attitude error can stay small. Rapid rotation can amplify sensitivities, and thrust limits can make a commanded correction infeasible. Those mechanisms need separate measurements. A threshold based only on angular speed or a Hessian sample does not establish whether the current task is controllable within the remaining time.
2.4 Geometric Interpretation
2.4.1 Moving Tangent Comparisons
A first variation is tangent to the family of model trajectories. For manifold-valued states, define errors in a local chart or with a specified retraction and compare tangent vectors consistently. On \(SO(3)\), for example, an attitude error may use the local logarithm of \(\bar R^TR\) with a declared branch. Its coordinates cannot be treated as a globally unique vector through arbitrary rotations.
The flow remainder depends on this representation. That dependence does not make it useless: a control design needs a metric and outputs. It does mean a small remainder in one set of coordinates is not automatically a coordinate-invariant statement about the motion’s difficulty.
2.4.2 Connection to Riemannian Curvature
Intrinsic Riemannian curvature requires a metric and its connection. A coordinate Hessian of a vector field is not that curvature tensor. Even motion on a flat one-dimensional line can have a nonlinear coordinate representation. On \(x>0\), take \(\dot x=1\) and introduce \(z=x^2\). Then
\[ \dot z=2\sqrt z,\qquad \frac{d^2\dot z}{dz^2}=-\frac1{2z^{3/2}}. \tag{14}\]
The state line remains intrinsically flat, yet this vector field has a nonzero coordinate Hessian. Its finite-time flow is
\[ z(t)=(\sqrt{z_0}+t)^2. \tag{15}\]
It is nonlinear in \(z_0\) despite the exact affine flow \(x(t)=x_0+t\). This counterexample is sufficient to reject a universal identification of flow remainder with intrinsic curvature. It also shows why the metric must be fixed before comparing “more nonlinear” trajectories.
2.4.3 Visual Example: Double Integrator and Pendulum
A double integrator with additive input is affine, so its exact flow has zero Taylor remainder for any admissible initial/input perturbation. Its numerical integration error depends on the numerical method. A pendulum has a nonlinear sine term, but its second derivative vanishes at the downward equilibrium. Around that equilibrium, the first omitted angular term is cubic. At other angles a quadratic term generally appears. Therefore “nonlinear system” does not imply a nonzero quadratic coefficient at every nominal point.
A direction reversal is a kinematic event where a chosen velocity component changes sign. It need not coincide with a large Hessian, a singularity, or maximum conversion of all system kinetic energy to potential energy. Other joints and the club can continue moving while one segment reverses.
3 Part II: Algorithms
4 Adaptive Timestep DDP
4.1 Motivation: Allocate Error Budgets Along a Motion
There are at least three different grids: the numerical integrator’s internal steps, the optimizer’s control/state mesh, and the physical controller’s update times. Refining one does not automatically refine the others. An embedded Runge–Kutta error estimate controls an integration defect. A Taylor remainder estimates a local model’s validity for a proposed perturbation. A controller update must also meet sensing, actuation, and computation deadlines.
An optimizer can use a large interval if the discrete flow over that interval is accurately integrated and differentiated and its local approximation is adequate. Conversely, a tiny interval does not guarantee a reliable first-order model for a large state/control update. Mesh density and trust-region size interact; neither replaces the other.
4.2 Adaptive Timestep Selection Rule
On one interval, suppose \(\|N\|\le q=M\rho^2/2\), where \(\rho\) bounds the actual joined deviation, and \(\|\Phi(t,s)\|\le e^{a(t-s)}\). A convenient Euclidean growth bound is any upper bound on the largest eigenvalue of \((A+A^T)/2\) throughout the interval. This follows by differentiating the squared norm of a variational solution. It accounts for non-normal transient growth that eigenvalues of \(A\) alone can miss.
For a comparison reset to \(r(0)=0\), the local error bound is
\[ R(h)=q\,\psi(a,h),\qquad \psi(a,h)=\begin{cases}(e^{ah}-1)/a,&a\ne0,\\h,&a=0.\end{cases} \]
For \(q>0\), budget \(\varepsilon_r\ge0\), and a non-negative growth bound \(a\), inversion gives
\[ h\le\begin{cases} \log(1+a\varepsilon_r/q)/a,&a>0,\\ \varepsilon_r/q,&a=0. \end{cases} \tag{16}\]
Use log1p and expm1 near zero. For \(q=0\), this particular local remainder condition imposes no interval limit. Other domain, integration, contact, and deadline conditions still apply. If a computed admissible step is below a hardware or optimizer minimum, report the unmet condition or change the plan; clipping upward does not preserve the bound.
The units check is decisive: \(q\) has units of inverse time, \(\varepsilon_r\) is dimensionless, and \(a\varepsilon_r/q\) is dimensionless. A square root of \(\varepsilon_r/q\) has units of square-root time and cannot be used as this interval rule.
With \(q=0.005\,\mathrm{s^{-1}}\), \(a=2\,\mathrm{s^{-1}}\), and \(\varepsilon_r=0.001\), the rule gives about \(0.168236\) s; without amplification it would give \(0.2\) s. These are manufactured comparison-equation parameters, not recommended golf or robot tolerances.
For multiple intervals, retain earlier error:
\[ R_{i+1}=e^{a_ih_i}R_i+q_i\psi(a_i,h_i),\qquad R_0=0. \]
Simply subdividing an unchanged comparison cannot erase its accumulated remainder. If relinearization resets a local comparison, the relationship to the original reference and its accumulated task error still needs to be tracked.
4.3 Algorithm 2.1: Refine a Fixed-Horizon Optimization Mesh
The following is an implementable design procedure, not a report of a validated new DDP solver:
- Specify terminal time \(T\), controls, state and input constraints, contact schedule, objective, integration tolerances, and task outputs. Start with positive intervals summing to \(T\).
- Roll out the nonlinear nominal trajectory using the same input interpolation and integration scheme whose discrete derivatives enter the optimizer. Include the quadrature rule in the stage cost.
- Solve the local optimization subproblem. Evaluate proposed state/control deviations, a regional enclosure when available, and propagated remainder indicators. Label sampled indicators as estimates.
- Accept or reject a candidate using its nonlinear rollout, actual cost decrease, constraint feasibility, and integration checks. A small remainder alone is insufficient. Limit update size when the local model predicts poorly.
- Split intervals that fail the chosen discretization or approximation criteria. Preserve total time exactly. Transfer controls using their specified interpolation, locate events separately, and reintegrate the nominal trajectory on the new mesh.
- Recompute derivatives, quadrature, terminal values, and constraints on that mesh. Reevaluate the old and new candidates with a common sufficiently accurate cost calculation before claiming improvement. Continue until the stated tolerances or budget terminate the procedure.
For a left-endpoint approximation to a continuous running cost \(\ell\), the stage term is \(h_i\ell(x_i,u_i)\); its state and control derivatives inherit \(h_i\). For a more accurate quadrature, differentiate that quadrature. Combining RK4 state rollout with unqualified Euler discrete Jacobians solves inconsistent local models.
4.4 Convergence Status and Assumptions
On a fixed mesh, iLQR approximates nonlinear dynamics to first order inside a second-order local cost model. Full DDP also includes contractions of the dynamics’ second derivatives with the value gradient. Both require an appropriate local solve, regularization, and accepted rollout; neither name guarantees convergence of arbitrary code. Li and Todorov’s biological movement paper develops iLQR and reports model-specific numerical comparisons; it is not the original DDP convergence proof (Li and Todorov 2004).
Changing the mesh changes the finite-dimensional optimization problem. Convergence of an inner solve does not establish convergence of the outer mesh sequence. A credible claim must specify consistent discretization, admissible mesh changes, termination, constraint treatment, regularity, and the norm in which solutions converge. ALTRO provides a concrete published approach to constrained trajectory optimization using augmented-Lagrangian and projection machinery; invoking its name does not validate the adaptive procedure above (Howell et al. 2019).
4.5 Residual-Triggered Mode Switching
4.5.1 Tracking LQR and Predictive Control
Tracking LQR can stabilize a nonlinear system locally around a nominal trajectory. MPC may use a linear or nonlinear prediction model. The relevant choice includes constraints, prediction horizon, model validity, available computation, and fallback feasibility; it is not simply “LQR for linear physics, MPC for nonlinear physics.”
4.5.2 Switching Logic
Choose an upper threshold \(\varepsilon_\uparrow\) and lower threshold \(\varepsilon_\downarrow<\varepsilon_\uparrow\) for a defined statistic. Count consecutive samples above the upper threshold and consecutive samples below the lower threshold. Reset the corresponding count whenever its condition fails; samples in the gap cannot count as consecutive evidence. Add a minimum dwell time when the design requires it.
Before transferring control, verify that the requested controller has a feasible current solution, compatible state/input conventions, and an admissible first action. Handle input slew limits and delayed solutions. If optimization misses its deadline, use an explicitly designed fallback within its valid region. A mode label is not a control policy.
Hysteresis reduces repeated threshold crossings but does not prove stability. A common Lyapunov function decreasing under every active controller can supply one route under its domain assumptions. Multiple Lyapunov functions require bounds on switching increases and sufficient decrease between switches. Constraints and contact transitions add their own conditions.
4.5.3 Computational Cost Analysis
Measure the whole loop: sensing, estimation, derivatives, optimization, integration, validation, transfer, and communication. Report hardware, solver tolerances, input size, warm starts, compilation policy, trial count, and worst-case or tail latency relevant to the deadline. Average optimizer time alone cannot establish real-time feasibility. No timing table in this article is presented as an observed speedup.
5 Tube MPC With Geometric Residuals
5.1 Residual-Based Tube Sizing
A tube encloses possible tracking errors under feedback and disturbances. A Taylor-validity neighborhood says where a local approximation is accurate. These are different sets. The first must be invariant or reachable under the stated dynamics; the second may merely impose an upper limit on perturbations.
5.2 Formulation
For nominal state \(\bar z(t)\) and error set \(E(t)\), write
\[ \mathcal T(t)=\bar z(t)\oplus E(t). \tag{17}\]
For a discrete local error model \(e^+=A_Ke+w\), a sufficient fixed-set condition is \(A_KE\oplus W\subseteq E\). A state constraint requires \(\bar z\in\mathcal X\ominus E\); under \(v=\bar v+Ke\), input admissibility requires \(\bar v\in\mathcal U\ominus KE\). The disturbance set must include the bounded nonlinear remainder and all other modeled uncertainty. A state-dependent remainder often makes finding \(E\) a self-consistent containment problem.
By contrast, in the nonexpansive, constant-bound case a local validity calculation gives only
\[ \rho\le\sqrt{\frac{2\varepsilon_r}{Mh}},\qquad M>0,\ h>0. \tag{18}\]
This square root is dimensionally valid as a deviation radius. It does not prove that the closed-loop error remains inside that radius.
5.3 Algorithm 2.2: Establish Containment Before Tightening Constraints
Consider the scalar teaching error dynamics \(\dot e=-2e+e^2+w\), \(|w|\le0.1\). For \(E=[-\rho,\rho]\), the worst outward velocity at the positive boundary is \(-2\rho+\rho^2+0.1\). At the negative boundary it is \(-2\rho-\rho^2+0.1\), so the first expression is the larger. A sufficient and, for this interval, exact boundary condition is
\[ \rho^2-2\rho+0.1\le0, \qquad 1-\sqrt{0.9}\le\rho\le1+\sqrt{0.9}. \]
Thus \(\rho=0.1\) gives inward boundary flow for every admissible disturbance, while \(\rho=0.01\) does not. Making a tube arbitrarily small fails because disturbances can push the state out. These are dimensionless manufactured dynamics and do not certify a humanoid or golfer.
For a practical construction: propose an error set; bound the nonlinear remainder over it; propagate the full uncertainty under the chosen feedback; verify containment and admissible inputs; then tighten nominal constraints and verify the initial error belongs to the set. If any check fails, revise the set, feedback, or nominal plan rather than asserting a tube from Taylor accuracy alone.
5.4 Tube-MPC Implication Under a Valid Bound
Mayne, Seron, and Raković study constrained linear discrete-time systems with bounded disturbances, including a nominal initial state chosen by the optimization and stability of a disturbance-invariant set (Mayne et al. 2005). This article’s scalar boundary calculation is an independent illustration of containment. Extending tube claims to nonlinear humanoid dynamics, changing contacts, and uncertain state estimates requires corresponding invariant/reachable-set and recursive-feasibility arguments; the linear result cannot be transferred by substituting a residual radius.
6 Part III: Applications
6.1 Quadrotor Aerobatics
6.1.1 Problem Setup
A rigid quadrotor has configuration \((p,R)\in\mathbb R^3\times SO(3)\) and velocities \((\dot p,\omega)\). The state manifold has dimension twelve; storing all nine entries of \(R\) produces eighteen numbers with six orthonormality constraints. A planar flip model has fewer states and cannot be reported as a full three-dimensional vehicle benchmark.
Specify mass, inertia, rotor allocation, thrust and rate limits, aerodynamic assumptions, and attitude-error convention. At a full flip, use consistent local attitude errors or manifold methods rather than a global small-angle Euler approximation.
6.1.2 Residual Analysis
The translational acceleration depends on thrust direction, so state-input cross terms matter. The rigid-body angular equation contains \(\omega\times I\omega\), whose sensitivities depend on inertia and rates. A large rotation angle alone neither gives a Hessian bound nor determines tracking difficulty. Compare predicted task errors in position and attitude, actuator saturation, and model mismatch under the same experiment.
6.1.3 Adaptive Timestep Evaluation
Use identical dynamics, initial/terminal conditions, constraints, cost, tolerances, and hardware for a fixed coarse mesh, fixed fine mesh, and adaptive mesh. Keep terminal time fixed. Reintegrate accepted solutions on a finer common reference before comparing terminal error and cost. Report all failed solves and deadline misses alongside successful trajectories. Separate derivative-compilation time from steady-state execution while reporting both.
6.1.4 Simulation Evidence
The checked calculations here concern the scalar flow and comparison equations. They establish the algebra used for residual budgeting, not a vehicle performance gain. A quadrotor evaluation needs a runnable model, solver configuration, initial conditions, saved trajectories, constraint residuals, and timing records. The evaluation above states what would make such a comparison meaningful; no synthetic percentage or unexplained table substitutes for those artifacts.
7 Humanoid Walking With Impacts
7.1 Problem Setup
Walking combines smooth constrained phases with contact events. In an ideal instantaneous impact model, positions remain continuous while generalized velocities can jump by a finite amount. A residual formed by comparing different contact modes at the same time can be large even for nearby initial states. This is a mismatch of hybrid comparisons, not an infinite physical state error.
For a fixed contact configuration, positive-definite mass matrix \(M_q\), independent contact-velocity constraints \(J_c\), no impulsive actuation, and a perfectly plastic sticking impact,
\[ \begin{aligned} M_q(\dot q^+-\dot q^-)&=J_c^T\Lambda,\qquad J_c\dot q^+=0,\\ \dot q^+&=\left[I-M_q^{-1}J_c^T (J_cM_q^{-1}J_c^T)^{-1}J_c\right]\dot q^-. \end{aligned} \tag{19}\]
The contact impulse must also satisfy the assumed unilateral and friction conditions. If it does not, the sticking reset is not admissible; slipping, restitution, or multiple-contact models require different equations. \(J_c\) denotes the full set of imposed velocity constraints, not automatically only the contact normal.
7.2 Residual Monitoring Across Contact Events
For a smooth guard \(g(t,z)=0\) and reset \(R(t,z)\), event-time variation contributes to the linearized jump. At a transverse, isolated event with the same nearby mode sequence, the saltation matrix is
\[ \Xi=D_zR+ \frac{(F^+-D_zR\,F^--D_tR)D_zg} {D_tg+D_zg\,F^-}. \]
Use smooth transition matrices within modes and \(\Xi\) at the event. The denominator must be nonzero. Grazing, simultaneous competing guards, and changes in event sequence require additional analysis. Kong and colleagues derive this event-time correction and its assumptions (Kong et al. 2024). A smooth Hessian bound should not be integrated through an instantaneous reset as though the vector field remained smooth.
For fall prediction, define the event to be predicted and a lead-time requirement, then compare calibrated prediction methods on held-out perturbations. Report false alarms, missed falls, lead times, contact-estimation uncertainty, and intervention outcomes. A residual spike after loss of balance is not evidence of useful early detection.
7.3 Experimental Validation Requirements
Hardware claims need an identified robot and controller, sensing and contact estimation, trial protocol, denominators, disturbances, failures, and data or traceable records. No physical ATRIAS experiment or fall-rate improvement is established here. Rigid-body teaching examples cannot establish muscle-level behavior or human fall prevention. Their value is in exposing which additional mechanisms and measurements a real evaluation needs.
7.4 Golf Swing Optimization
7.4.1 Problem Setup
For golf, define the output before choosing a residual norm. A useful impact-state output can include clubhead velocity, face orientation, path direction, and strike location relative to the ball. Ball launch speed, direction, and spin require a collision model and ball/club parameters; maximizing a club-speed coordinate alone does not determine them.
The state model may include joint configurations and velocities, club bending/torsion, activation states, and tendon or soft-tissue states when those mechanisms are relevant. Their inclusion changes both the drift and input maps. Joint torques are not direct observations of individual muscle forces, and a torque-input model does not by itself explain neural commands or activation delay. Validate each model level against the measurements that can identify it.
7.4.2 Residual-Aware Trajectory Refinement
Let \(y=c(z(T))\) be a smooth pre-impact task output and \(C_T=Dc(\bar z(T))\). To first order, the contribution of the flow remainder to the output discrepancy is \(C_Tr(T)\). The norm contribution generated at time \(s\) is weighted by \(\|C_T\Phi(T,s)\|\), not merely by \(\|N(s)\|\):
\[ \|C_Tr(T)\|\le \int_{t_0}^T\|C_T\Phi(T,s)\|\,\|N(s)\|\,ds. \]
For a fully rigorous output bound, also enclose the nonlinear remainder of \(c\) over the relevant states. For a differentiable event-to-output map, include event-time and reset sensitivities. If the impact sequence changes, a single smooth map may cease to describe the comparison.
This explains how apparently separate swing questions intertwine. Proximal motion changes the configuration and velocities from which later distal motion develops. Those changes alter moment arms, inertial coupling, control authority, and the sensitivity of impact outputs. Club compliance adds stored deformation and delayed recovery; activation dynamics restrict which torques can be changed before impact. The eventual consequence depends on this evolving system, the available time, and the output direction—not on a universal rule that all proximal motion must slow or that every local remainder must be minimized.
One can make this distinction precise with an adjoint. For a scalar terminal score \(j(z(T))\), solve \(-\dot p=A^Tp\) backward with \(p(T)=\nabla j\). The leading signed contribution of a small unmodeled forcing \(d(t)\) is \(\int p(t)^Td(t)\,dt\). Magnitude-only bounds protect against unknown directions; signed sensitivity explains a particular modeled perturbation. Cancellation in a signed integral must not be promoted to a worst-case guarantee. For a hybrid motion, the adjoint also crosses the appropriate transpose jump maps.
A sound refinement study therefore varies numerical mesh, perturbation size, model complexity, and task metric separately. It checks whether a claimed improvement survives a finer integration, a richer club/body model, plausible parameter uncertainty, and repeated initial conditions. It reports both successful and failed trials. A phase-based mesh must sum to the stated total swing duration; the short collision phase should be treated with a resolved contact law or a justified impact map.
An expert golfer might permit variation that scarcely changes the shot while controlling other directions tightly. That is a testable task-level hypothesis. It does not establish that the nervous system computes Hessians, minimizes a coordinate residual, or implements an LQR/MPC switch. Learning, perception, fatigue, noise, and unobserved physiological states need their own evidence.
8 Part IV: Implementation
9 Derivatives and Reproducible Examples
Automatic differentiation can supply derivatives of a specified smooth program. It does not certify that the program models the golfer, that a sampled derivative bounds a neighborhood, or that a nonsmooth contact branch has an ordinary Hessian. For discrete optimal control, differentiate the actual discrete rollout and cost; for a continuous-time Taylor bound, differentiate the continuous field in the declared coordinates.
9.1 Reproducible Bound Calculations
The repository implements the comparison recurrence in src/affine_control/residuals.py and the two scalar teaching functions in src/tools/residual_control_examples.py. Run this compact example from the repository root with its Python dependencies installed:
import numpy as np
from src.affine_control.residuals import predict_residual_bound
from src.tools.residual_control_examples import residual_step_limit, scalar_nonlinear_flow
initial, horizon = 0.01, 1.0
terminal = scalar_nonlinear_flow(initial, np.array([horizon]))[0]
actual_remainder = terminal - initial * np.exp(horizon)
old_expression = initial**2 * np.expm1(2 * horizon) / 2
conditional_bound = predict_residual_bound(
np.array([2.0]), np.array([terminal]), np.array([horizon]), np.array([1.0])
)
assert old_expression < actual_remainder < conditional_bound
step = residual_step_limit(forcing=0.005, budget=0.001, growth=2.0)predict_residual_bound takes one interval upper bound per entry for the Hessian and actual joined deviation, non-negative interval durations, and propagator growth rates. Omitting growth rates selects zero growth, so the legacy sum requires a nonexpansive propagator. Negative growth rates are supported in propagation; the interval-inversion helper uses the simpler non-negative-growth convention. The calculations use ordinary floating point, not verified interval arithmetic.
The historical compute_hessian_bound name denotes a finite-difference point estimate of the largest component state-Hessian norm, with control held fixed. It neither includes mixed derivatives nor certifies a regional vector bound. Its corrected documentation states this explicitly. ResidualMonitor reports a scaled measurement/reference discrepancy and a hysteretic label. It implements neither an observer nor a feasible MPC controller.
9.2 Performance Comparison
The regression suite checks the comparison formula against independent quadrature, preserves accumulation under subdivision, verifies the exact scalar flow against numerical ODE integration, and tests the interval rule by integrating its comparison equation. It also checks the pendulum Hessian sign and the scalar tube’s boundary condition. Run:
python -m pytest tests/test_residual_flow_bounds.py tests/test_residual_control_examples.pyThese tests validate teaching calculations and API contracts. They do not measure optimizer speed, physical tracking accuracy, or clinical outcomes. The repository’s explicitly labeled DDP mock cannot provide evidence of those outcomes. A publishable adaptive-solver result requires the complete solver and evaluation artifacts described above.
9.3 Practical Guidelines
9.3.1 When to Use Residual-Aware Methods
Use the analysis when a local approximation influences a consequential decision and you can define its error, propagation, and domain. It can guide offline mesh refinement, perturbation limits, and model diagnostics. A fixed-rate controller may still use residual diagnostics without changing its physical update frequency. Hybrid models can use the analysis within smooth modes and event sensitivities at transitions.
9.3.2 Choosing a Residual Tolerance
Allocate the acceptable task error among estimation, model discrepancy, linearization, integration, and actuation. Propagate each contribution to the same task coordinates before comparing them. For correlated random errors, include cross-covariances; adding variances requires the appropriate independence assumptions. For deterministic bounded errors, a triangle-inequality sum is conservative but clear. A noise floor does not automatically set an acceptable impact error, and a dimensionless tolerance has no portable meaning without its scales.
9.3.3 Debugging Adaptive Algorithms
| Observation | Check Before Changing the Algorithm |
|---|---|
| Large predicted bound but small observed model remainder | Conservatism of the regional and propagator bounds; metric scaling |
| Small estimate but large nonlinear-model discrepancy | Domain coverage, actual versus linearized deviation, missing mixed terms |
| Accurate model integration but large measurement innovation | Estimation, omitted forces, parameters, timestamps, sensor calibration |
| Cost improves after regridding | Common horizon, input interpolation, quadrature, feasibility, reference integration |
| Large discrepancy near contact | Event time, active constraints, reset admissibility, nearby mode sequence |
| Frequent controller switches | Counter resets, hysteresis, dwell time, uncertainty calibration, feasible transfer |
| Required step below the allowed minimum | An unmet budget; changing the minimum cannot prove the original bound |
10 Conclusion
11 Summary
The useful object is a chain of reasoning: a defined model comparison generates a remainder; the dynamics propagate it; the output map determines its task consequence; and a controller acts only within its feasibility and stability conditions. In golf and humanoid motion, that chain passes through inertial coupling, actuation, compliance, estimation, and contact. Treating any one residual as the whole explanation breaks the chain. Keeping those relationships explicit makes local linear analysis a strong component of nonlinear modeling.
11.1 Connections to Other Articles
Use the series index to connect first variations and integration, DDP/iLQR, contraction, and hybrid dynamics. The complete golf-swing chapter develops the mechanical context. The present argument supplies an error-budgeting layer connecting those models to their task outputs.
12 Further Reading
- Li and Todorov (2004): iLQR for nonlinear biological movement models, with model-specific comparisons (Li and Todorov 2004).
- Howell, Jackson, and Manchester (2019): ALTRO’s constrained trajectory-optimization formulation (Howell et al. 2019).
- Mayne, Seron, and Raković (2005): robust MPC for constrained linear systems with bounded disturbances (Mayne et al. 2005).
- Kong and colleagues (2024): saltation matrices and the assumptions behind hybrid first variations (Kong et al. 2024).
A simplified model predicts a swing only approximately. This article asks where the prediction misses, how that miss grows or shrinks, and whether it changes the shot. For a longer plain-language walk-through, read the accessible explanation.
A Miss Is a Clue
The gap between prediction and result tells us to look closer. It does not name the cause. A timing slip, a wrong force model, or a noisy sensor can all look alike.
Errors Travel Forward
An error made early is carried through the rest of the motion. It can grow or fade. A small wrist error can matter a lot at the clubface, while a bigger error elsewhere may barely change the shot.
What Is Not Shown
The article works out error bounds and worked examples in simple cases. It does not show that adaptive methods are faster, prevent falls, or help golfers. Those remain proposed tests.
13 Appendix: Mathematical Proofs
13.1 Proof of Theorem 1.1
Fix time \(t\) and define \(\varphi(s)=F(t,\bar z+se,\bar v+s\delta v)\). Applying the fundamental theorem of calculus twice gives
\[ \varphi(1)-\varphi(0)-\varphi'(0) =\int_0^1\int_0^s\varphi''(\sigma)\,d\sigma\,ds =\int_0^1(1-\sigma)\varphi''(\sigma)\,d\sigma. \]
With \(w=(e,\delta v)\), \(\varphi''(\sigma)=D^2F[w,w]\) along the segment. The bilinear norm bound and \(\int_0^1(1-\sigma)d\sigma=1/2\) establish \(\|N\|\le M\|w\|^2/2\). This argument needs \(C^2\) regularity on the segment. A pointwise second-order Taylor polynomial has an \(o(\|w\|^2)\) remainder under \(C^2\); a uniform \(O(\|w\|^3)\) bound requires stronger regularity, such as a locally Lipschitz Hessian.
Since \(\dot e=Ae+B\delta v+N\) and \(\dot\eta=A\eta+B\delta v\), subtraction gives \(\dot r=Ar+N\). Multiplying by a fundamental-matrix inverse, integrating, and multiplying back yields \(r(t)=\Phi(t,t_0)r(t_0)+\int\Phi(t,s)N(s)ds\). Here \(r(t_0)=0\) by definition. Taking norms proves Equation 6. There is no step in which the propagator may be replaced by one without an additional bound.
For a perturbation family with \(\|e(t_0)\|+\|\delta v\|_\infty\le\alpha\), bounded first derivatives on a common finite-horizon domain give \(\|e\|_\infty\le C\alpha\) by a Grönwall comparison. A uniformly bounded Hessian and propagator then give \(\|r\|_\infty\le C'\alpha^2\). This supplies the assumptions behind the quadratic order claim and shows why its constants depend on the horizon and domain.
For closed-loop feedback \(\delta v=Ke\), the open-loop comparison above has \(\dot\eta=A\eta+BKe\), using the actual input difference. A different comparison uses \(\dot\eta_K=(A+BK)\eta_K\). Define \(r_K=e-\eta_K\); then \(\dot r_K=(A+BK)r_K+N(t,e,Ke)\). Its propagator is the closed-loop one. For a nonlinear policy, perform Taylor expansion of the composed closed-loop field, including all policy derivatives. The two comparisons are individually valid but their symbols and constants cannot be interchanged.
13.2 Derivation of the Adaptive Interval Rule
Assume \(\|N\|\le q\) and \(\|\Phi(t,s)\|\le e^{a(t-s)}\) within an interval of duration \(h\). With zero initial remainder,
\[ \|r(h)\|\le q\int_0^he^{a(h-s)}ds=q\psi(a,h). \]
For \(a>0\) and \(q>0\), \(q(e^{ah}-1)/a\le\varepsilon_r\) is equivalent to \(e^{ah}\le1+a\varepsilon_r/q\), giving Equation 16. The continuous limit as \(a\to0\) is \(h\le\varepsilon_r/q\). For \(a<0\), the bound increases toward \(q/|a|\). If \(\varepsilon_r\ge q/|a|\), this condition allows any finite \(h\) within the validity domain; otherwise the logarithmic inverse still applies with a positive logarithm argument. The helper deliberately requires \(a\ge0\) to keep this case distinction out of its public contract.
For nonzero incoming bound \(R_i\), variation of constants gives the recurrence in the main text. Iterating it propagates every interval’s contribution through all later intervals. Holding the same coefficients and forcing while subdividing gives the same exact comparison solution. Consequently, local mesh refinement is useful only through changes in approximation, integration, control representation, or verified domains; subdivision alone cannot improve an unchanged continuous-time bound.
13.3 Exercises and Checks
- For \(F(z,v)=zv\), expand about \((\bar z,\bar v)\) and identify the mixed remainder. Explain why \(D^2_{vv}F=0\) does not remove it.
- Derive the scalar counterexample’s exact solution by separation of variables. Find its pole and check the small-\(\epsilon\) ratio to the unpropagated expression.
- For \(A=\begin{bmatrix}-1&10\\0&-1\end{bmatrix}\), compute the transition matrix. Explain why two negative eigenvalues do not imply \(\|\Phi(t,0)\|_2\le1\).
- Verify the interval rule by solving \(\dot R=2R+0.005\), \(R(0)=0\), at the stated budget. Convert the entire calculation from seconds to milliseconds.
- Derive the two outward boundary velocities of the scalar robust-tube example. Explain both the lower and upper admissible radii.
- A golf model predicts clubface orientation at a fixed time, while measured ball contact occurs earlier. Identify the additional event sensitivity needed before attributing the mismatch to swing control.
Answer Checks: (1) \(N=e\delta v\). (2) The pole is \(\log(1+1/\epsilon)\); the limiting ratio is \(2e^T/(e^T+1)\). (3) \(\Phi=e^{-t}\begin{bmatrix}1&10t\\0&1\end{bmatrix}\), which can transiently amplify norms. (4) \(h\approx0.168236\) s, or \(168.236\) ms with both rates divided by \(1000\). (5) The positive boundary is limiting and gives \(1\pm\sqrt{0.9}\); too little room cannot contain disturbances, while a sufficiently large positive error overwhelms the stabilizing linear term. (6) Account for the contact guard’s time shift, the corresponding state/output change, and any reset before comparing equivalent events.
References
Citation
@online{olson2026,
author = {Olson, Dieter},
title = {Residual-Aware {Control:} {Bounds} and {Decisions} for
{Nonlinear} {Motion}},
date = {2026-01-18},
url = {https://affinedrift.com/articles/tangent-hyperplane-articles/Advanced/Residual-Aware_Control.html},
langid = {en}
}
Computational Overhead and Benefit
Derivative estimates and repeated optimization may cost more than they save, especially under tight deadlines.
Biomechanics and Model Assumptions
Rigid-body residuals cannot explain every difference between a golfer and a simulated mechanism.
The Danger of Trusting Residuals
A large innovation may come from noise or a missing force rather than local nonlinearity.