Forward Dynamics, Counterfactuals, and Drift
A moving club carries the consequences of everything that brought it to its present state. Removing a torque does not remove that history. Position, velocity, shaft deformation, muscle activation and contact conditions determine what happens next, together with the inputs that remain. A counterfactual simulation asks a precise question about changing one part of that continuation. It can show how a specified intervention changes an acceleration, a trajectory, an energy balance or an impact condition. Those are different measurements, and their distinctions are essential to understanding a golf swing.
The useful design question is how to establish a state from which the desired outcome requires manageable, feasible corrections. Passive dynamics can help carry a movement toward that outcome. They can also carry it away. This chapter develops the calculations needed to distinguish these possibilities, using a completely specified two-link example before returning to golf and humanoid modeling. No population-level claim about muscular effort follows from the algebra alone.
Defining the Intervention
Control-Affine Dynamics and the Meaning of Zero
For a declared state \(x\) and input \(u\), write
\[ \dot x=F(x,u)=f(x)+G(x)u. \]
Here \(f(x)=F(x,0)\) is the drift and \(G(x)u\) is the same-state input contribution. At a fixed state, the dependence on \(u\) is affine. Both contributions obey physics; the terminology does not partition the motion into voluntary and inevitable parts. If the physical model contains explicit time dependence, write \(f(x,t)\) and \(G(x,t)\), or augment the state with a clock. We use autonomous notation where that is adequate.
The identity
\[ \frac{\partial f}{\partial u}=0 \]
states a property of this representation. It is not evidence that observational data identify the physical effect of removing an actuator. Indeed, the input change \(u=\alpha(x)+\beta(x)w\), with invertible \(\beta\), gives drift \(f+G\alpha\) and input matrix \(G\beta\). Zero \(w\) is then the physical policy \(u=\alpha(x)\), not zero \(u\). A control coordinate is part of the intervention definition.
In an ideal torque-actuated rigid-body model, setting a joint motor torque to zero is well defined. In a muscle-driven model, zero neural excitation, zero activation and zero muscle force are different interventions. Existing activation can evolve after excitation is removed; elastic tissue and the shaft can retain stored energy. A prescribed grip trajectory can continue to supply boundary power even when some internal joint torques are set to zero. State which channel is removed and which dynamics, supports and inputs remain. The word “passive” should only be used with that declaration.
Four Counterfactual Objects
The pointwise zero-torque counterfactual (ZTCF) evaluates \(f(x(t))\) at the actual state. A stitched ZTCF collects such evaluations along the actual trajectory. It is a diagnostic field sampled on that trajectory, not generally a trajectory of \(f\).
A forward ZTCF starts from one specified state and integrates zero declared input:
\[ \dot x^{(0)}=f(x^{(0)}),\qquad x^{(0)}(t_b)=x(t_b). \] A branched ZTCF repeats this forward experiment at different branch times \(t_b\). Every branch inherits the actual history only up to its own start. Its subsequent state must be solved independently. Plotting several branch results at a common horizon answers a different question from restarting the branch at every plotted sample.
These definitions permit exact same-state input recovery:
\[ \dot x(t)-f(x(t))=G(x(t))u(t). \]
But with \(\Delta x=x-x^{(0)}\), a forward branch obeys
\[ \Delta\dot x=f(x)-f(x^{(0)})+G(x)u. \]
The drift changes as the two states separate. Integrating the same-state input contribution alone therefore does not recover their separation. Removing several actuators in different orders can give different marginal trajectory effects; a sum of marginal attributions requires an explicit convention for distributing those interactions.
For smooth mechanical dynamics \(\dot q=v\), \(\dot v=a_0(q,v)+B_a(q,v)u\), both branches have the same position and velocity at \(t_b\). For a right-continuous constant input near that time,
\[ \begin{aligned} \Delta v(t_b+h)&=B_a(x_b)u_bh+O(h^2),\\ \Delta q(t_b+h)&=\tfrac12 B_a(x_b)u_bh^2+O(h^3). \end{aligned} \]
A small short-horizon position gap is partly a consequence of this second-order onset. It does not establish negligible acceleration, velocity change or impact error. A motor can influence a sensitive output while barely changing a coarse plot of the whole trajectory.
When a Counterfactual Exists
Local Solutions and Their Domain
On an open domain, a locally Lipschitz drift has a unique local solution through each initial state. A continuously differentiable drift satisfies that local regularity requirement. Its solution extends to a maximal interval while it remains in the domain and does not escape every compact subset. Smoothness alone does not guarantee an arbitrarily long simulation: \(\dot x=x^2\) has solution \(x(t)=x_0/(1-x_0t)\) and blows up at \(1/x_0\) for \(x_0>0\).
For \(f\in C^k\), \(k\geq1\), the local flow is \(C^k\) in the initial state; a trajectory is \(C^{k+1}\) in time. Continuity alone is insufficient for uniqueness. Contact switching, frictional transitions, actuator saturation and impacts may require piecewise-smooth or differential-inclusion models instead of this smooth theorem.
Flow Maps and Linearized Propagation
Write \(\phi_\tau(x_b)\) for the autonomous zero-input flow over elapsed time \(\tau\). Then
\[ \phi_0=\mathrm{id},\qquad \phi_{\tau+s}(x)=\phi_\tau(\phi_s(x)) \]
where all expressions are defined. Negative-time inversion is also local; a global group needs a complete vector field. Prescribed time-dependent inputs instead require a two-time evolution map.
The state-transition matrix (STM) of the drift solves
\[ \dot\Phi_0(t,t_b)=D f(x^{(0)}(t))\Phi_0(t,t_b),\qquad \Phi_0(t_b,t_b)=I. \]
It gives first-order initial-state sensitivity, not the exact difference between finite perturbations. Generally it is a time-ordered product of local propagators, not the exponential of an averaged Jacobian. Eigenvalues at one instant do not determine its singular values or finite-time amplification. This distinction becomes decisive when an intervention changes both the configuration and the remaining time to impact.
Mechanical Drift, Velocity Resets, and Constraints
A Consistent Force Convention
For autonomous natural mechanics with independent coordinates, use
\[ \begin{aligned} M(q)\dot v+C(q,v)v+g_q(q)+d(q,v)&=B(q)u,\\ g_q:=\nabla_q V,\qquad a_0&=-M^{-1}(Cv+g_q+d),\\ B_a&=M^{-1}B. \end{aligned} \]
Assume \(M\) is positive-definite on the coordinate domain. The holding-force term \(g_q\) is on the left; gravity’s applied generalized force is \(-g_q\). If \(V\) includes springs, \(g_q\) includes elastic forces as well. The dissipative term \(d\) is on the left and has \(v^Td\geq0\) for a passive damper. Forces excluded from this equation, such as a moving boundary’s reaction or prescribed aerodynamic load, must be introduced before calling it a model of a particular swing.
At a fixed state, multiplying by \(M^{-1}\) is linear. For the stated force partition this yields
\[ a_0=a_{\mathrm{grav}}+a_{\mathrm{cor}}+a_{\mathrm{damp}},\qquad a_{\mathrm{grav}}=-M^{-1}g_q,\quad a_{\mathrm{cor}}=-M^{-1}Cv,\quad a_{\mathrm{damp}}=-M^{-1}d. \] Here “grav” means gravity alone only when \(V\) contains only gravitational potential. This is an acceleration decomposition in declared coordinates, not a unique partition of energy transferred between bodies. Individual entries of \(C\) are coordinate dependent, and \(C\) itself need not be the unique matrix factor of its bias vector \(Cv\).
What a Zero-Velocity Counterfactual Removes
A pointwise zero-velocity counterfactual (ZVCF) replaces the velocity by zero at fixed configuration and specified internal states. In the simple model above, with \(C(q,0)0=0\) and \(d(q,0)=0\), the zero-input version is
\[ a_{\mathrm{ZVCF}}=a_0(q,0)=-M(q)^{-1}g_q(q). \]
A control-preserved velocity reset instead gives \(a(q,0,u)=a_0(q,0)+B_a(q,0)u\). Label which version is used. “Zero velocity” is not “static equilibrium”: this acceleration can be nonzero.
In this natural mechanical model,
\[ a_0(q,v)-a_0(q,0)=-M^{-1}\{C(q,v)v+d(q,v)-d(q,0)\}. \]
Thus the difference includes damping whenever it is present. More general velocity-dependent forces, actuator characteristics and internal states add their own terms. A reset may also change kinetic energy discontinuously; it is a diagnostic intervention, not a free physical operation.
Contact Reactions Also Respond to Input
Suppose a fixed smooth contact mode imposes \(A\dot v=-b\) and the force equation is \(M\dot v=Bu-h+A^T\lambda\), with \(h=Cv+g_q+d\). For independent constraint rows, define
\[ \begin{aligned} \Lambda&=(AM^{-1}A^T)^{-1},\\ P&=I-M^{-1}A^T\Lambda A,\\ a_0^c&=-PM^{-1}h-M^{-1}A^T\Lambda b,\\ B_a^c&=PM^{-1}B. \end{aligned} \]
These expressions follow by solving for \(\lambda=-\Lambda\{b+AM^{-1}(Bu-h)\}\). A contact reaction therefore generally changes with \(u\); freezing the actual reaction while deleting torque defines a different, usually inconsistent experiment.
For a moving holonomic constraint \(\varphi(q,t)=0\), admissible velocities satisfy \(Av+\varphi_t=0\). Setting \(v=0\) can violate this relation. Specify an admissible reset or limit the ZVCF to a formal off-manifold diagnostic. Unilateral contact additionally requires valid force signs and friction conditions; if a branch lifts off or slips, switch models and check the new mode. A fixed-base reduced model can have no explicit multipliers while still representing physical constraints. A golfer’s ground and grip conditions cannot be inferred from the appearance of the reduced equations.
Finite-Horizon and Impact Counterfactuals
Inputs Act Through the Remaining Trajectory
Along a smooth nominal trajectory with prescribed input \(u(t)\), let
\[ A(t)=D f(x(t))+\sum_j u_j(t)D G_j(x(t)). \]
With no initial-state perturbation, the first-order response to a small input variation is
\[ \delta x(t_f)=\int_{t_b}^{t_f}\Phi(t_f,s)G(x(s))\delta u(s)\,ds, \] where \(\dot\Phi=A\Phi\). For feedback \(u=\pi(x,t)\) plus an added perturbation, include \(G D_x\pi\) in \(A\). A replayed command and a feedback controller that adapts to the branch are different experiments. A failed actuator can also change friction or electrical loading; its model must say whether it is free, braked or otherwise altered.
For output \(y=h(x(t_f))\), premultiply the integral by \(D h(x(t_f))\). A pulse of amplitude \(\delta u_b\) and duration \(\Delta\) has leading response \(D h\Phi(t_f,t_b)G(x_b)\delta u_b\Delta\). Changing a bounded input at one isolated time has zero effect on an ordinary ODE solution; an impulse is a different model. Singular values of a pulse-response matrix require input and output weights. For distributed inputs, the finite-horizon controllability Gramian measures an \(L^2\)-energy problem, not automatically peak-bounded physiological feasibility.
One rigorous small-gap statement uses a Lipschitz bound \(L\) for \(f\) on a region containing both branches and a uniform bound \(\|G(x(t))u(t)\|\leq\epsilon\). Gronwall’s inequality gives
\[ \|\Delta x(t_b+h)\|\leq \begin{cases} \epsilon(e^{Lh}-1)/L,&L>0,\\ \epsilon h,&L=0. \end{cases} \]
These are absolute bounds in a declared norm. A ratio with a large drift numerator supplies neither a small \(\epsilon\) nor a small propagation factor. For a smooth output with Lipschitz constant \(L_h\), multiply the state-gap bound by \(L_h\). A finite removal of a large torque generally requires nonlinear branch integration; the linear sensitivity is an approximation to be checked by perturbation refinement.
Impact Is an Event, Not Merely a Clock Time
Let impact be the transversal event \(\psi(x(t_e),t_e)=0\). If \(\delta x^-\) is the fixed-time pre-event perturbation, differentiation of the event equation yields
\[ \delta t_e=-\frac{D_x\psi\,\delta x^-}{\psi_t+D_x\psi F^-},\qquad \delta x_e=\delta x^-+F^-\delta t_e. \]
For \(y=h(x_e,t_e)\),
\[ \delta y=D_xh\,\delta x^-+(D_xh F^-+h_t)\delta t_e. \]
A post-impact comparison also differentiates the reset law and propagates its resulting perturbation. The denominator must be nonzero; near grazing, event time can be highly sensitive, and a perturbation can change whether contact happens at all.
For the elementary motion \(q(t)=vt\) toward \(q=1\), the fixed-time derivative \(\partial q/\partial v=t\) is nonzero. At the crossing event the position derivative is zero, while \(dt_e/dv=-1/v^2\). This is why comparing clubhead positions at the nominal impact time can answer a different question from comparing clubface orientation and speed at each branch’s own ball encounter. Both are legitimate when labeled.
A Reproducible Attribution Protocol
First define the output and tolerance: for example, pre-impact clubhead velocity and orientation at an explicit contact event. Specify the model, coordinates, internal states, contact law, input channels and feasible limits. Reconstruct the initial state and nominal input with uncertainty, checking that the forward model reproduces the measured motion over the relevant interval.
Next run the declared branches. Record branch time, horizon, remaining command or feedback policy, solver tolerances and any mode changes. Evaluate pointwise acceleration contributions separately from forward trajectory gaps. Compare task errors, event times, actuator work and constraint residuals. Refine step size or tolerances and perturb uncertain parameters; check whether the conclusion survives those changes.
Finally state the supported claim at its actual scope. A small modeled impact change under one torque deletion is evidence about that deletion from that state. It is not a measurement of muscle inactivity, the source of earlier energy, or how a different swing should be coached. An infeasible branch or a missed ball encounter is an outcome to report, not a reason to silently compare a different event.
Energy and Parameter Interventions
Energy Is a Balance Over a Boundary
For the autonomous natural model with fixed ideal supports,
\[ E=T+V=\tfrac12v^TM(q)v+V(q),\qquad \dot E=u^TB^Tv-\mathcal D,\quad \mathcal D=v^Td\geq0. \] The identity \(v^TCv=\tfrac12v^T\dot Mv\) cancels the coordinate-dependent inertial terms in total power. Motor torque is power-conjugate to \(B^Tv\), which is not necessarily the vector of absolute segment angular rates. Moving supports, prescribed external loads and explicit time dependence require their corresponding power terms. At a stationary, no-slip ground contact, contact force can change momentum while the contact point has zero mechanical power.
For zero input and no dissipation, \(E^{(0)}\) is constant. With passive damping it is nonincreasing, not strictly decreasing at every instant. Conservation does not imply periodic oscillation: free translation, rotation and separatrix motion are counterexamples. Nor does conservation imply that all nearby trajectories contract.
For actual and zero-input trajectories starting from the same state,
\[ E(t)-E^{(0)}(t)=\int_{t_b}^t u^TB^Tv\,ds -\int_{t_b}^t(\mathcal D-\mathcal D^{(0)})\,ds. \]
Add the difference of boundary-power terms if present. The two branches usually dissipate different amounts because their states differ. Even in the conservative case, the gap measures signed net actuator work, not metabolic cost or the sum of positive work. Positive and negative actuator work can cancel while control remains consequential.
To infer energy received by an individual link, draw its boundary and evaluate interface wrench power using that body’s interface velocity. Equal and opposite wrenches cancel in a whole-system balance only when their conjugate velocities agree, with relative motion and deformation accounted for. Off-diagonal mass-matrix entries alone do not identify which link receives the most energy. Prescribing torso deceleration also requires whatever forces and power enforce that prescription; the torso does not become a free source by being treated as an input trajectory.
Changing Gravity, Mass, or an Actuator
Let a parameter \(p\) enter the complete vector field \(F(x,u,p)\). With fixed prescribed input, the parameter sensitivity \(S=\partial x/\partial p\) satisfies
\[ \dot S=F_xS+F_p,\qquad S(t_b)=\partial_p x_b. \] If the command depends explicitly on \(p\), add \(F_u u_p\); if it is feedback, use the closed-loop state derivative as well. Differentiating only the instantaneous inverse mass matrix misses the propagated state change.
Changing gravity at a fixed configuration in a fixed mode can leave \(B_a\) unchanged while changing \(a_0\). It does not imply a linear superposition of the new trajectory and the old one. Changing mass or geometry can change \(M,C,V,B\) and contact behavior together. Such parameter dependence does not destroy affinity in the input, but it changes the counterfactual model and requires recomputation. A local sensitivity is useful for small perturbations; a substantial equipment change needs a new finite simulation and a check of whether the same event sequence persists.
Identifying Drift From Data
The Separation Must Be Identifiable
When the analytical model is uncertain, a chosen basis can represent \(f\) and \(G\) in a regression
\[ \dot x_k=\Theta(x_k,u_k)\theta+r_k. \]
For a model constrained to be control-affine, include state features and state features multiplied by each input; do not insert arbitrary quadratic input terms and still call the result affine. A unique noiseless fit of unrestricted linear coefficients requires a full-column-rank stacked regressor. Practical estimation additionally needs adequate conditioning, representative states and a defensible model for residuals and measurement errors.
Feedback data can conceal the distinction. If observations use only \(u=\pi(x)\), then for any compatible matrix field \(H(x)\),
\[ \widetilde f=f+H\pi,\qquad \widetilde G=G-H \]
produce the same observed vector field. Yet their zero-input predictions differ. Repeated or nearby states under varied inputs help distinguish the effects, and structural physical constraints can reduce the ambiguity. Exact state repetition is not necessary for every parametric model; assess the actual regression rank and experiment design. Persistent excitation is meaningful relative to the chosen regressors and the estimation theorem, not a generic assurance that an input happened to vary.
Noise, Validation, and Model Choice
For independent position noise of variance \(\sigma^2\), a forward-difference velocity estimate has noise variance \(2\sigma^2/\Delta t^2\). Faster sampling alone can increase derivative noise. Smoothing, integral formulations or explicit measurement models can help, but their bias and time resolution must be assessed. Inverse-dynamics torques and differentiated kinematics can share errors, so an ordinary regression may suffer errors-in-variables and feedback correlation.
Sparse identification with control fits a selected library of state and input features; Kaiser, Kutz and Brunton provide a concrete formulation and numerical demonstrations (Kaiser et al. 2018). Its usefulness depends on the library and data, and its reported model comparisons do not establish universal superiority. Linear ARX, state-space realization and DMD with control have different representation assumptions. Balanced truncation reduces an existing suitable model; it does not by itself identify a nonlinear biomechanical drift.
Validate on withheld trajectories and input policies, including relevant zero-input or safely approximated interventions where evidence is available. Check energy and constraint residuals, report uncertainty, and distinguish interpolation from extrapolation into states absent from training. A compelling actual-trajectory fit is necessary for some uses but cannot alone validate an unobserved counterfactual. Treat the resulting golf interpretation as a model prediction until independent evidence supports that intervention.
Stability and Optimal Control Are Separate Questions
A Contraction Statement With a Complete Proof
Let \(f\) be continuously differentiable and \(J_f=D f\). On a convex forward-invariant domain, suppose solutions exist for all times of interest and
\[ \tfrac12(J_f(x)+J_f(x)^T)\preceq-\alpha I,\qquad \alpha>0 \]
uniformly throughout that domain. Then any two solutions there contract in Euclidean distance. To prove this, put \(\delta=x_1-x_2\) and use the exact mean-value identity
\[ \begin{aligned} f(x_1)-f(x_2)&=\int_0^1J_f(x_2+s\delta)\delta\,ds,\\ \frac{d}{dt}\tfrac12\|\delta\|^2&\leq-\alpha\|\delta\|^2,\\ \|\delta(t)\|&\leq e^{-\alpha(t-t_b)}\|\delta(t_b)\|. \end{aligned} \] Convexity keeps the connecting segment in the domain; forward invariance keeps both solutions there. No discarded Taylor remainder is needed. Pointwise negative eigenvalues without a uniform margin are insufficient: \(\dot x=-x^3\) on \(x>0\) has negative derivative but only algebraic decay toward zero.
For a state-dependent metric \(P(x,t)\), the differential condition is \(\dot P+J_f^TP+PJ_f\preceq-2\alpha P\), with the derivative evaluated along the drift. Uniform positive lower and upper bounds on \(P\) yield uniform equivalence to Euclidean distance. Finite-distance conclusions require admissible connecting paths and flow-domain conditions. The metric formulation and its path argument are developed by Lohmiller and Slotine (Lohmiller and Slotine 1998); the kinetic-energy mass matrix is not automatically a contraction metric.
Eigenvalues, Hamiltonian Motion, and Feedback Margins
The constant matrix \(J=\left[\begin{smallmatrix}-1&10\\0&-1\end{smallmatrix}\right]\) has two negative eigenvalues, yet \(e^{Jt}=e^{-t}\left[\begin{smallmatrix}1&10t\\0&1\end{smallmatrix}\right]\) amplifies some initial errors before they decay. Its symmetric part has largest eigenvalue \(4\). For time-varying systems even asymptotic stability cannot be inferred from instantaneous eigenvalues alone.
Conservative canonical Hamiltonian flow preserves phase-space volume. It cannot uniformly contract every full-state direction for all time on a regular region in a uniformly equivalent metric. Some directions can contract transiently while others expand. Damping, feedback, projection onto an output, or transverse stability of an orbit are different properties; specify which one is being measured. In an undamped swing model, an observed narrowing of clubhead paths is not by itself a full-state contraction certificate.
Small feedback can preserve an existing stability margin under explicit conditions. For an equilibrium \(f(x^*)=0\) and a perturbation \(r(x)=G(x)K(x-x^*)\), a Euclidean contraction bound remains valid if the symmetric part of \(D r\) has norm below \(\alpha\) throughout the relevant invariant domain. The target remains an equilibrium because \(r(x^*)=0\). Local exponential stability also survives sufficiently small smooth perturbations of a Hurwitz equilibrium Jacobian. Neither statement says that an arbitrary target is an equilibrium or that all small gains accomplish a specified finite-time task.
A Large Drift Ratio Does Not Determine Optimal Effort
For an unconstrained input and a smooth finite-horizon value function \(\mathcal V(t,x)\) with running cost \(\ell(x)+u^TRu\), \(R\succ0\), the HJB minimization gives
\[ \begin{aligned} 0&=\mathcal V_t+\min_u\{\ell(x)+u^TRu+\nabla_x\mathcal V^T(f+Gu)\},\\ u^*&=-\tfrac12R^{-1}G^T\nabla_x\mathcal V. \end{aligned} \] Appropriate terminal data and solution conditions remain necessary. With feasible input bounds, minimize over that set; for a convex set this is an \(R\)-weighted projection of the unconstrained minimizer. A large \(\|f\|\) says nothing by itself about the value gradient or whether the drift points toward a low-cost outcome.
An explicit counterexample comes from \(\dot x=ax+u\) and \(J=\int_0^\infty(x^2+u^2)\,dt\). Write \(\mathcal V=Px^2\). Substitution gives
\[ P^2-2aP-1=0,\qquad P=a+\sqrt{a^2+1},\qquad u^*=-Px. \] For \(a=10\), \(P\approx20.05\): strong unstable drift calls for strong stabilizing feedback. For \(a=-10\), \(P\approx0.04988\): favorable stable drift reduces the optimal gain for this particular objective. The sign, horizon, input cost and target matter. The Riccati/HJB construction is standard (Tedrake 2024); these scalar examples expose the missing implication directly.
This unconstrained LQR example does not assert feasibility under a fixed torque bound. Indeed, with \(a>0\) and \(|u|\leq u_{\max}\), any \(x>u_{\max}/a\) has \(\dot x>0\) even under maximum opposing input and cannot be regulated to zero. Large drift relative to capacity can mean loss of recoverability. Even stable drift can require sustained input at a nonzero target: \(\dot x=-x+u\) needs \(u=1\) to hold \(x=1\). Exploiting drift is therefore an optimization and feasibility question, not a universal zero-control theorem.
Worked Example: An Underactuated Double Pendulum
Geometry, Inertia, and Torque Coordinates
Two massless rigid rods of lengths \(l_1,l_2\) carry point masses \(m_1\) at the first endpoint and \(m_2\) at the second endpoint. The base pivot is fixed; motion is planar. The absolute angles \(q_1,q_2\) are measured from the downward vertical in the same positive sense. This is a point-mass model, distinct from a model with uniformly distributed rod mass.
Let \(\delta=q_1-q_2\), \(k=m_2l_1l_2\). Endpoint velocities give
\[ M(q)=\begin{bmatrix}(m_1+m_2)l_1^2&k\cos\delta\\k\cos\delta&m_2l_2^2\end{bmatrix},\qquad V=-(m_1+m_2)gl_1\cos q_1-m_2gl_2\cos q_2. \]
The determinant is \(m_2l_1^2l_2^2(m_1+m_2\sin^2\delta)>0\) for positive masses and lengths. The Christoffel factor and holding-force vector are
\[ C=\begin{bmatrix}0&k\sin\delta\,v_2\\-k\sin\delta\,v_1&0\end{bmatrix},\qquad g_q=\begin{bmatrix}(m_1+m_2)gl_1\sin q_1\\m_2gl_2\sin q_2\end{bmatrix}. \] In particular, the off-diagonal Coriolis signs are opposite. Direct differentiation verifies \(\dot M-2C\) is skew-symmetric. At a small positive angle with zero rates, gravity initially accelerates a simple hanging pendulum back toward zero, consistent with the positive \(g_q\) on the equation’s left.
Here a single motor acts between the base and the first rod, with torque \(u=\tau_1\); the second joint is unactuated, so \(B=(1,0)^T\). If an internal elbow motor were added with relative rotation \(q_2-q_1\), its generalized torque column would be \((-1,1)^T\). The full-state input vector for the single base motor is the four-by-one column \((0,0,(M^{-1})_{11},(M^{-1})_{21})^T\), not a stack of two rows of \(M^{-1}\).
For \(m_1=m_2=1\,\mathrm{kg}\), \(l_1=l_2=1\,\mathrm m\) and \(g=9.81\,\mathrm{m/s^2}\),
\[ \begin{aligned} M&=\begin{bmatrix}2&\cos\delta\\\cos\delta&1\end{bmatrix},& M^{-1}&=\frac1{1+\sin^2\delta}\begin{bmatrix}1&-\cos\delta\\-\cos\delta&2\end{bmatrix},\\ g_q&=\begin{bmatrix}19.62\sin q_1\\9.81\sin q_2\end{bmatrix},& G&=\frac1{1+\sin^2\delta}\begin{bmatrix}0\\0\\1\\-\cos\delta\end{bmatrix}. \end{aligned} \] The displayed numerical mass entries have units \(\mathrm{kg\,m^2}\); holding torques are in \(\mathrm{N\,m}\) and the acceleration entries of \(G\) convert torque into angular acceleration. The state equations are \(\dot q=v\), \(\dot v=M^{-1}(B\tau_1-Cv-g_q)\).
Two Integrated Branches and Their Actual Evidence
Start both branches at \(q_1=\pi/4\), \(q_2=\pi/6\), \(v_1=2\,\mathrm{rad/s}\) and \(v_2=-1\,\mathrm{rad/s}\). Apply constant \(\tau_1=5\,\mathrm{N\,m}\) to the actual branch and zero torque to the ZTCF for two seconds. There is no damping or impact. Angles below are continuously unwrapped in this chart.
| Time (s) | \(q_1\) | \(q_2\) | \(v_1\) | \(v_2\) | \(a_{0,1}\) | \(5(M^{-1})_{11}\) | \(\rho_1\) |
|---|---|---|---|---|---|---|---|
| 0.0 | 0.7854 | 0.5236 | 2.0000 | -1.0000 | -9.742 | 4.686 | 2.079 |
| 0.5 | 0.9514 | 0.3060 | -1.4024 | 0.1886 | -10.707 | 3.672 | 2.916 |
| 1.0 | -0.0516 | 0.2870 | -0.8654 | -2.2171 | 4.953 | 4.503 | 1.100 |
| 1.5 | -0.0242 | -1.0481 | -0.1023 | -1.0901 | -2.871 | 2.891 | 0.993 |
| 2.0 | -0.0517 | -0.5357 | -0.3469 | 3.1249 | -6.586 | 4.110 | 1.602 |
Actual Trajectory With Constant Base Torque. Angles Are in Radians, Rates in Radians per Second, and Accelerations in Radians per Second Squared. The Final Column Is the First-Coordinate Ratio.
| Time (s) | \(q_1\) | \(q_2\) | \(v_1\) | \(v_2\) | \(\|\Delta q\|\) |
|---|---|---|---|---|---|
| 0.0 | 0.7854 | 0.5236 | 2.0000 | -1.0000 | 0.0000 |
| 0.5 | 0.5263 | 0.6641 | -2.7217 | 1.2024 | 0.5558 |
| 1.0 | -0.2760 | -0.2033 | 0.4717 | -4.6111 | 0.5393 |
| 1.5 | -0.4161 | -1.3513 | -0.5783 | 0.2386 | 0.4955 |
| 2.0 | -0.5652 | -0.2543 | 0.6854 | 2.8132 | 0.5855 |
Zero-Torque Trajectory From the Same Initial State. The Final Column Is the Norm of Its Position Difference From the Actual Branch, in Radians.
The actual table’s \(\rho_1=|a_{0,1}|/|5(M^{-1})_{11}|\) is a first-coordinate realized-input ratio. It also equals that component’s capacity ratio if the single motor limit is exactly \(|\tau_1|\leq5\,\mathrm{N\,m}\). It is not the full-vector \(\rho_W\). It rises from \(2.079\) to \(2.916\), then falls to \(1.100\) and \(0.993\) before reaching \(1.602\). There is no monotonic progression to negligible control. The position gap is \(0.5558\) rad at \(0.5\) s and \(0.5855\) rad at \(2\) s; these are coordinate-norm separations, not a fraction of motion attributable to an actuator.
At \(1\) s on the actual branch, the acceleration decomposition is
\[ a_{\mathrm{grav}}=\begin{bmatrix}3.271\\-5.862\end{bmatrix},\qquad a_{\mathrm{cor}}=\begin{bmatrix}1.682\\-1.835\end{bmatrix},\qquad a_0=\begin{bmatrix}4.953\\-7.698\end{bmatrix} \quad\mathrm{rad/s^2}. \]
The base torque contributes approximately \(4.503\,\mathrm{rad/s^2}\) to the first acceleration. It acts through both acceleration components, with coupling that changes sign or magnitude with configuration. The unactuated coordinate is not disconnected from the input.
The ZTCF conserves \(E=T+V\). The actual branch obeys the independent identity
\[ E(t)-E(0)=5\{q_1(t)-q_1(0)\}. \]
At two seconds this is approximately \(-4.1854\) J. Positive torque has done negative net work because the net first-angle displacement is negative. The two branches consequently do not share an energy shell. Their separation need not accumulate monotonically, and this short trajectory provides no global bound on unwrapped angle differences.
Reproducing and Checking the Numbers
The tabulated values are generated by the repository’s fixed-step RK4 solver using \(5{,}000\) steps over two seconds; its Coriolis matrix is derived from finite differences of the specified mass matrix. Independent validation uses the analytic bias \(C v=(\sin\delta\,v_2^2,-\sin\delta\,v_1^2)^T\) with adaptive DOP853 integration. Agreement at the printed precision, tolerance refinement and the work-energy identity test different possible errors. Agreement between two solvers is useful evidence only when they implement the same declared physical model, and it does not validate a human swing model.
Run this example from the repository root with Python 3.12 and the project dependencies.
import numpy as np
from scipy.integrate import solve_ivp
from scripts.generate_worked_examples import (
ch07_mass_matrix, ch07_grad_potential,
)
from src.affine_control.dynamics import christoffel_coriolis
initial = np.array([np.pi / 4, np.pi / 6, 2., -1.])
times = np.linspace(0, 2, 5)
for torque in (0., 5.):
def rate(t: float, state: np.ndarray) -> np.ndarray:
q, v = state[:2], state[2:]
bias = christoffel_coriolis(ch07_mass_matrix, q, v) @ v
rhs = np.array([torque, 0.]) - bias - ch07_grad_potential(q)
return np.r_[v, np.linalg.solve(ch07_mass_matrix(q), rhs)]
for tolerance in (1e-9, 1e-11):
result = solve_ivp(
rate, (0, 2), initial, method="DOP853",
t_eval=times, rtol=tolerance, atol=tolerance / 100,
)
if not result.success:
raise RuntimeError(result.message)
q1, q2, v1, v2 = result.y
potential = -9.81 * (2 * np.cos(q1) + np.cos(q2))
energy = (v1**2 + v2**2 / 2
+ np.cos(q1 - q2) * v1 * v2 + potential)
work = torque * (q1 - q1[0])
print(torque, tolerance, np.max(abs(energy - energy[0] - work)))
print(result.y.T)To reproduce a finite torque removal at a later time, first integrate the actual state to that time and initialize a new zero-torque solve there. Do not substitute the original zero-input branch’s state, which carries a different history. To compare impact outcomes, add the intended contact geometry and event detection before interpreting the result as a striking motion.
How the Pieces Fit Together in a Golf Swing
Preparation, Transport, and Correction
Earlier muscle and ground interactions establish momentum, configuration and stored elastic energy. Configuration determines inertia, moment arms, coupling and which actuator directions can change the task. Evolving velocity changes both kinematic transport and some generalized bias forces. These facts explain why a late snapshot cannot identify how the present motion was produced.
Counterfactual branching measures what additional inputs do from that snapshot onward. A short remaining interval can limit changes in clubhead position, while orientation, speed or event timing remain sensitive in particular directions. Existing motion can therefore carry much of the path while selected inputs still matter for impact. Neither observation contradicts the other: they concern different outputs and time scales.
At the top of the backswing, the club’s reversal does not require every body segment, muscle state and elastic mode to have zero velocity. Gravitational acceleration depends on the actual configuration and supports, and need not resemble the subsequent controlled motion. During the fast phase, a large inertial or radial force does not by itself show that the useful tangential or orientational control direction is weak. The size of a wrist torque is also not the same quantity as its work, a grip force, or a muscle’s metabolic demand.
A Testable Version of the Effort-to-Ease Idea
A precise hypothesis is that a particular preparation policy places the system in states from which the desired impact can be achieved with lower subsequent effort while satisfying error, actuator and contact constraints. To test it, compare policies at a declared impact objective and evaluate their feasible continuations. Measure the remaining work or chosen effort cost, sensitivity to uncertainty, and reserve for correcting perturbations. A low nominal effort solution with no correction reserve may be fragile.
Use several interventions to expose the connections: remove a specified torque after different branch times; retain or alter the feedback of the other actuators; perturb segment inertia within measured uncertainty; and compare the resulting event-conditioned task errors. These experiments can distinguish a favorable inherited state from a favorable instantaneous drift direction, and both from active correction. An apparent advantage that disappears when grip motion is no longer prescribed may have depended on uncounted boundary actuation.
For humanoid motions, the same logic connects balance, contact feasibility and task execution. Ground contact can redirect momentum without positive work at an ideal stationary contact, yet the internal actuation needed to maintain or exploit that contact can be substantial. For golf, any claim about a typical athlete or coaching intervention requires measurements and validation beyond this two-link illustration. The model can sharpen the experimental question and show which quantities must be measured.
Chapter Synthesis
Same-state decomposition isolates a declared input’s instantaneous contribution. A forward branch reveals the consequence of removing it over a specified interval. An energy balance accounts for work and storage across a defined boundary. Directional and finite-horizon sensitivities connect inputs to a task, and event sensitivities account for changes in when contact occurs. Stability and optimal control require their own metric, objective and feasibility assumptions.
Together these tools support a constructive account of skilled movement: establish useful states, understand how geometry and dynamics transport them, and retain enough feasible authority to meet the task under uncertainty. The strength of the argument comes from connecting those calculations and testing their assumptions, rather than treating a single ratio as a complete explanation of the swing.
Exercises
- A Nonlinear Oscillator and Its Branch. For the nondimensional model \(\ddot z+z^3=u\), start at \((z,\dot z)=(1,0)\). Integrate zero input and the feedback \(u=-2z\). Plot positions, velocities and \(E=\dot z^2/2+z^4/4\). Verify \(\dot E=u\dot z\) and conservation of \(E+z^2\) on the feedback branch. With declared capacity \(|u|\leq3\), compute \(\rho=|z|^3/3\) along each branch, separately from the realized-input ratio. Explain why a zero actual input makes the latter undefined or infinite without implying actuator failure.
- Velocity Resets and Constraints. For Exercise 1 compute zero-input and control-preserved ZVCF accelerations at the initial state; compare them with the actual acceleration. There is no gravity in this nondimensional oscillator: identify the restoring term. Next impose the moving constraint \(q-t=0\) on a unit mass. Show why setting its admissible velocity to zero violates the constraint, and explain what additional reset rule a physical counterfactual would require.
- Conditional Quadratic Scaling. At fixed configuration, write the double-pendulum Coriolis acceleration for \(v=sw\). State conditions on \(q,w\) and actuator capacity under which \(\rho_W/s^2\) has a positive finite limit. Give a direction or configuration where the quadratic coefficient vanishes. Compare with the constant-length single pendulum and explain the role of Cartesian radial acceleration.
- Work, Potential, and Dissipation. For fixed ideal supports, derive \(\Delta T=W_u-\Delta V-W_d\), where \(W_u=\int u^TB^Tv\,dt\) and \(W_d=\int\mathcal D\,dt\). Derive the actual-minus-ZTCF total-energy balance. Verify the two-link constant-torque identity numerically. Construct a nonzero control history with zero net work and explain why it need not leave the trajectory unchanged.
- Identifying an Input Channel. Use the model \(\ddot q=\theta_1\sin q+\theta_2\dot q+\theta_3u\). Generate noiseless and noisy data under two distinct input policies; examine the rank and singular values of the regressor. Fit parameters on one collection and test on the other. Construct a policy that makes a pair of columns dependent. Repeat with differentiated noisy positions and explain any bias or instability in the recovered parameters.
- Stability in a Declared Metric. Compute the undamped double-pendulum Jacobian near its hanging equilibrium and compare the full-state flow with a damped linearization. Explain why canonical conservative flow cannot uniformly contract every state direction in a uniformly equivalent metric. For \(J=[-1,10;0,-1]\), calculate transient amplification and find a constant positive-definite \(P\) satisfying \(J^TP+PJ\prec0\).
- A Torque Pulse and an Impact Event. For a pendulum with a small constant torque pulse over \([t_b,t_b+\Delta]\), integrate the variational response to a fixed terminal time and compare it with finite differences of nonlinear trajectories as pulse amplitude decreases. Repeat for a transversal angle-crossing event, including its time shift. State the input/output weights before identifying the most influential pulse direction. Explain why a change at a single isolated time is different from a pulse or an impulse.
- A Three-Link Energy Attribution. Specify a torso–arm–forearm model, joint coordinates, applied torques, boundary constraints and segment inertias. If torso deceleration is prescribed, solve for the required enforcing reaction. Compute each link’s energy change and interface wrench power, including relative joint motion. Compare a torque-removal branch with the actual motion. Explain which results follow from off-diagonal inertia terms and which require the complete force, velocity and boundary-power calculation.
References
Citation
@online{untitled,
author = {},
url = {https://affinedrift.com/articles/The_Geometry_of_Motion/quarto/ch07_counterfactuals.html},
langid = {en}
}