Applications: Biomechanics, Robotics, and Beyond
A Mechanically Consistent Golf Teaching Model
Coordinates, States, and Rigid-Body Inertia
Consider a planar serial chain with three revolute coordinates \(q=(q_1,q_2,q_3)\). The first angle is measured from the horizontal; the second and third are relative joint angles. The absolute link orientations are \(\theta_i=\sum_{j=1}^i q_j\). Three configuration degrees of freedom require six first-order states \((q,\dot q)\). Adding two bending amplitudes \(\eta=(\eta_1,\eta_2)\) requires ten states \((q,\eta,\dot q,\dot\eta)\), not a six-state system with a different mass matrix.
Let \(l_i\) be link lengths, \(r_i\) the distances from each proximal joint to its center of mass, \(m_i\) the masses, and \(I_i\) the moments of inertia about those centers. With translational and angular Jacobians \(J_{v,i}\) and \(J_{\omega,i}\), \[ T_r=\frac12\sum_{i=1}^3\left(m_i\|J_{v,i}\dot q\|^2+I_i(J_{\omega,i}\dot q)^2\right),\qquad M_r=\sum_{i=1}^3\left(m_iJ_{v,i}^{\mathsf T}J_{v,i}+I_iJ_{\omega,i}^{\mathsf T}J_{\omega,i}\right). \] This expression includes both translation and rotation. Replacing each link by a point mass changes the model and must be stated. For independent coordinates and positive link inertias, nonzero joint velocity produces positive kinetic energy; a singular end-point position Jacobian does not make this joint-space inertia indefinite.
An explicit formula provides a useful independent check. Define \[ \begin{aligned} a&=I_1+I_2+I_3+m_1r_1^2+m_2(l_1^2+r_2^2)+m_3(l_1^2+l_2^2+r_3^2),\\ b&=I_2+I_3+m_2r_2^2+m_3(l_2^2+r_3^2),\qquad d=I_3+m_3r_3^2,\\ h&=l_1(m_2r_2+m_3l_2)\cos q_2,\\ j&=m_3l_1r_3\cos(q_2+q_3),\qquad k=m_3l_2r_3\cos q_3. \end{aligned} \] Then \[ M_r(q)=\begin{bmatrix} a+2(h+j+k)&b+h+j+2k&d+j+k\\ b+h+j+2k&b+2k&d+k\\ d+j+k&d+k&d \end{bmatrix}. \] The cross terms are not optional corrections: they express the velocity of a distal mass produced by several joints simultaneously. In particular, the coefficient of \(k\) in the first two entries includes contributions from both participating velocities. Energy computed from this matrix must agree with energy computed directly from the moving centers of mass.
Use \(m=(10,5,0.81)\) kg, \(l=(0.30,0.35,1.15)\) m, \(r_i=l_i/2\), and \(I=(0.25,0.08,0.089)\) kg m\(^2\). These represent a rigid carrier for the calculation; labeling its coordinates as particular anatomical joints does not validate that identification. At \(q=(60,-30,-20)\) degrees, the reproducible result is \[ M_r=\begin{bmatrix}2.7750 & 1.3863 & 0.5998 \\ 1.3863 & 0.9955 & 0.5100 \\ 0.5998 & 0.5100 & 0.3568\end{bmatrix}\ \mathrm{kg\,m^2}. \] Its determinant is approximately \(0.0682\) in the corresponding cubed inertia units; its smallest eigenvalue is approximately \(0.0501\) kg m\(^2\). These are checks on this numerical model, not clinical or performance measures.
A Flexible Appendage Must Contribute All Its Kinetic Energy
For an illustrative flexible extension, add a beam of mass \(m_s=0.30\) kg along the third link. This mass is additional to the specified rigid carrier. It must not also be counted inside \(m_3\) when interpreting the combined model. No separate clubhead mass moving with the flexible tip is included here. A physical club would require that mass and its rotational inertia, along with an appropriate shaft model and grip boundary condition.
Let \(s\in[0,l_3]\), \(t=(\cos\theta_3,\sin\theta_3)\), and \(n=(-\sin\theta_3,\cos\theta_3)\). Write a beam material point as \[ p(s,q,\eta)=p_h(q)+s\,t(q)+n(q)\sum_{a=1}^2\phi_a(s)\eta_a. \] At the undeformed configuration \(\eta=0\), its velocity has the form \(J_q(s)\dot q+J_\eta(s)\dot\eta\), with \(J_\eta=n[\phi_1\ \phi_2]\). The same material-point integral must generate every block: \[ \begin{aligned} A&=M_r+\int_0^{l_3}\rho_l J_q^{\mathsf T}J_q\,ds,\\ C&=\int_0^{l_3}\rho_l J_q^{\mathsf T}J_\eta\,ds,\qquad D=\int_0^{l_3}\rho_l J_\eta^{\mathsf T}J_\eta\,ds,\\ M_f&=\begin{bmatrix}A&C\\C^{\mathsf T}&D\end{bmatrix},\qquad \rho_l=m_s/l_3. \end{aligned} \] Including \(C\) and \(D\) while omitting the beam contribution to \(A\) does not describe that kinetic energy. It can even produce an indefinite matrix for otherwise positive physical parameters. Such a failure is an assembly error, not evidence of a kinematic singularity.
With dimensionless mode shapes and amplitudes in meters, \(A\) has units kg m\(^2\), \(C\) has units kg m, and \(D\) has units kg. Positive definiteness is preserved by invertible coordinate changes, but raw eigenvalues and condition numbers of the full mixed-unit matrix depend on the amplitude scales. Compare them only after declaring those scales.
For the displayed posture, numerical integration gives \[ A=\begin{bmatrix}3.2056 & 1.7293 & 0.8220 \\ 1.7293 & 1.2780 & 0.6990 \\ 0.8220 & 0.6990 & 0.4891\end{bmatrix},\qquad C=\begin{bmatrix}0.1594 & -0.0496 \\ 0.1368 & -0.0371 \\ 0.0981 & -0.0157\end{bmatrix},\qquad D\simeq\begin{bmatrix}0.0750 & 0.0000 \\ 0.0000 & 0.0750\end{bmatrix}. \] The units follow the preceding definitions. Rounding makes \(D\) appear exactly diagonal; the underlying quadrature has a small off-diagonal residual.
The shape functions are the first two clamped-free, uniform Euler–Bernoulli beam shapes without a tip mass. They can be derived from \(EI\phi''''=\rho_l\omega^2\phi\) with \(\phi(0)=\phi'(0)=\phi''(l_3)=\phi'''(l_3)=0\). With \(z=\beta_a s/l_3\), \[ \begin{aligned} \widetilde\phi_a(s)&=\cosh z-\cos z-\sigma_a(\sinh z-\sin z),\\ \sigma_a&=\frac{\cosh\beta_a+\cos\beta_a}{\sinh\beta_a+\sin\beta_a},\qquad \cosh\beta_a\cos\beta_a=-1,\\ \phi_a(s)&=\widetilde\phi_a(s)/\widetilde\phi_a(l_3). \end{aligned} \] The first two positive roots are approximately \(1.875104\) and \(4.694091\). Dividing by the signed tip value makes both \(\phi_a(l_3)=1\); dividing by a common positive constant leaves the second mode with the opposite sign. Either sign convention is physically acceptable if the amplitude and every coupling term transform consistently. It is inconsistent to claim a unit positive tip while using the other convention.
Here the reduced spring matrix is defined separately as \(K_\eta=\operatorname{diag}(D_{11}\,40^2,D_{22}\,120^2)\), approximately \(\operatorname{diag}(120,1080)\) N/m. It is symmetric and positive. These are prescribed reduced-model frequencies in rad/s. They do not correspond to one uniform beam stiffness: that beam would have \(\omega_2/\omega_1=(\beta_2/\beta_1)^2\simeq6.267\), rather than \(3\). A measured nonuniform shaft, head inertia, axial loading, and rotation may require different modes and additional terms. The displayed inertia and springs alone are not a complete flexible forward-dynamics simulator.
Eliminate Accelerations Without Losing the Bias Terms
For a fully specified flexible model, collect gravity, velocity-dependent forces, elastic forces, and damping in \(h_q\) and \(h_\eta\), using a consistent sign convention: \[ A\ddot q+C\ddot\eta+h_q=Bu,\qquad C^{\mathsf T}\ddot q+D\ddot\eta+h_\eta=0. \] Here \(C\) denotes an inertia cross block. It is distinct from the conventional Coriolis matrix in the product \(\mathbf{C}\dot q_{\text{sys}}\), where \(q_{\mathrm{sys}}=(q,\eta)\). A notation change cannot remove that velocity-dependent force from the model.
Define \(S_q=A-CD^{-1}C^{\mathsf T}\), \(H=S_q^{-1}\), and \(\Gamma=-D^{-1}C^{\mathsf T}\). Direct substitution gives \[ \begin{aligned} \ddot q&=H(Bu-h_q+CD^{-1}h_\eta),\\ \ddot\eta&=-D^{-1}h_\eta+\Gamma\ddot q. \end{aligned} \] Thus the acceleration input blocks are \(HB\) and \(\Gamma HB\). The definition of \(\Gamma\) contains no \(H\); including \(H\) there and multiplying by it again counts the mobility twice. The drift also contains the modal bias \(h_\eta\). Removing it would remove elastic restoring forces and any other unactuated modal forces assigned to that term.
For the ten-state ordering \(x=(q,\eta,\dot q,\dot\eta)\), define \(a_0=H(-h_q+CD^{-1}h_\eta)\). The first-order affine model is \[ f(x)=\begin{bmatrix}\dot q\\\dot\eta\\a_0\\-D^{-1}h_\eta+\Gamma a_0\end{bmatrix},\qquad G(x)=\begin{bmatrix}0\\0\\HB\\\Gamma HB\end{bmatrix}. \] Each row has the appropriate derivative units. In particular, angles cannot be subtracted from gravity torques inside a state derivative.
At the example posture, \[ S_q=\begin{bmatrix}2.8340 & 1.4142 & 0.6032 \\ 1.4142 & 1.0103 & 0.5123 \\ 0.6032 & 0.5123 & 0.3574\end{bmatrix},\qquad H=\begin{bmatrix}1.3860 & -2.7604 & 1.6179 \\ -2.7604 & 9.1216 & -8.4168 \\ 1.6179 & -8.4168 & 12.1323\end{bmatrix}. \] Positive definiteness follows structurally: for any nonzero \(v\), substitute the velocity pair \((v,-D^{-1}C^{\mathsf T}v)\) into the positive kinetic-energy quadratic form to obtain \(v^{\mathsf T}S_qv>0\). This proof requires the full, consistently assembled \(M_f\succ0\). Numerical inversion cannot rescue an inconsistent kinetic energy.
A Rigid Counterfactual Is a Different Model
The supplied zero-torque trajectory uses the six-state rigid chain, with \[ M_r(q)\ddot q+c_r(q,\dot q)+g_r(q)=u,\qquad a_{\mathrm{drift}}=-M_r^{-1}(c_r+g_r). \] It uses \(M_r\), not the flexible Schur complement. Mixing the latter with rigid-only bias terms would generally break the energy balance. With conservative gravity and \(u=0\), the rigid model conserves \(E=\tfrac12\dot q^{\mathsf T}M_r\dot q+V_r(q)\). Its fixed base can exert reaction forces even when the generalized joint torques are zero.
Initialize \(q=(60,-30,-20)\) degrees and \(\dot q=(12,-18,-25)\) rad/s. The absolute angular velocities are \((12,-6,-31)\) rad/s; interpreting the three relative rates as independent measured segment speeds would be incorrect. At this state, \[ \begin{aligned} g_r&=(30.25,\; 14.34,\; 4.50)^{\mathsf T}\ \mathrm{N\,m},\\ c_r&=(120.25,\; 11.14,\; -17.42)^{\mathsf T}\ \mathrm{N\,m},\\ a_{\mathrm{drift}}&=(-118.58,\; 70.47,\; 134.83)^{\mathsf T}\ \mathrm{rad/s^2}. \end{aligned} \] The velocity-force contribution exceeds gravity in the first component at this state. That observation neither establishes dominance in every coordinate nor proves a universal transition during human downswing.
A fourth-order Runge–Kutta integration over \(0.10\) s uses \(400\) steps of \(0.00025\) s. Selected outputs are shown below. “Rigid tip” is the endpoint of this chain, not a separately modeled flexible clubhead.
| \(t\) (s) | \(q_3\) (rad) | \(\dot q_3\) (rad/s) | Rigid-tip speed (m/s) | \(\|a_{\mathrm{drift}}\|_2\) (rad/s²) |
|---|---|---|---|---|
| 0.000 | -0.3491 | -25.00 | 35.37 | 192.9 |
| 0.050 | -1.2301 | -8.29 | 33.26 | 473.8 |
| 0.100 | -0.9462 | 24.21 | 21.66 | 900.9 |
The tip speed decreases from \(35.37\) to \(21.66\) m/s in this particular calculation while the acceleration norm increases. Those two facts can coexist because acceleration changes direction as well as speed, and because energy redistributes between kinetic and potential terms. This table supplies no controlled comparison and no evidence that a golfer should relax their muscles. Numerical credibility requires step refinement and an energy check; physiological credibility would additionally require validated masses, contacts, passive tissues, muscle forces, and initial conditions.
From Drift Decomposition to a Testable Golf Argument
Declare the Intervention and the Quantity Being Compared
At a fixed state of a torque-driven rigid model, acceleration separates exactly as \(\ddot q=a_0(q,\dot q)+M_r^{-1}u\). This is additivity in the input at one state. It does not imply that the trajectories produced by two different input histories add, because their evolving states change both terms.
A dimensionless drift-control ratio can be defined for a chosen acceleration output: \[ \rho_a=\frac{\|W_a a_0\|_2}{\|W_a M_r^{-1}u\|_2}. \] Here \(W_a\) specifies component scales or a declared projection and normalization. If the denominator is zero and the numerator is nonzero, the ratio is infinite; if both vanish, it is undefined. A numerical epsilon changes the definition and must be reported. The unforced trajectory above therefore cannot be assigned a finite drift-control ratio from its acceleration norm alone.
The full first-order drift \(f=(\dot q,a_0)\) combines rates and accelerations. Its Euclidean norm is not \(\|a_0\|\), and its comparison with \((0,M_r^{-1}u)\) depends on arbitrary units unless a state metric is specified. A metric-based ratio of full vector fields is a different quantity. Even acceleration itself changes nontrivially under nonlinear coordinates: for \(z=\psi(q)\), \[ \ddot z=D\psi(q)\ddot q+D^2\psi(q)[\dot q,\dot q]. \] The second term belongs to the transformed zero-input acceleration. Consequently, a large coordinate ratio is not an invariant amount of physiological effort or mechanical energy.
Velocity-dependent forces are often quadratic in velocity at fixed posture. Under \(\dot q\mapsto s\dot q\), the rigid \(c_r\) scales as \(s^2\), but it can vanish or cancel gravity in particular components. Meanwhile the applied torque and posture can change with speed. No theorem therefore makes this ratio increase through every swing or makes a particular numerical threshold a universal phase boundary.
A finite counterfactual should instead specify the same initial state, model parameters, contact rules, and future input intervention, then compare a stated output. “Set generalized active torques to zero” is one mathematical intervention. Removing neural feedback, reducing activation, removing all muscle force, and removing contact forces are different interventions. Already developed muscle force and an ongoing motor command do not disappear merely because a new sensory correction cannot arrive before impact. A meaningful claim about late control must include the actuator and delay model.
Acceleration Sequences Are Not Energy-Transfer Measurements
For the rigid chain, actuator power is \(u^{\mathsf T}\dot q\) and the total energy balance is \(\dot E=u^{\mathsf T}\dot q\) under the conservative assumptions above. For a single segment, power also enters through joint forces acting at moving points and moments acting with that segment’s absolute angular velocity. At an ideal joint with a shared point velocity, equal and opposite joint-force powers cancel in the combined energy balance while transferring energy between the individual segments. If the application-point velocities differ, as in a sliding contact, their sum need not vanish. Joint moments have distinct proximal and distal powers whose sum corresponds to the relative-coordinate actuator power.
A proximal angular-speed decrease can accompany distal acceleration, but it does not, by itself, identify where the distal kinetic energy came from. Translational energy, potential energy, active work, and energy exchange with other segments can all matter. To establish transfer, define the bodies, reference points, sign conventions, and interval; calculate the force and moment powers; integrate them; and check each segment’s energy balance. The temporal order of speed peaks is useful descriptive information, not a substitute for that accounting.
This distinction also resolves an apparent paradox in the teaching trajectory. A large acceleration norm can coexist with decreasing endpoint speed, while mechanical energy stays constant. A velocity component, an endpoint speed, a segment kinetic energy, and the total system energy are four different observables. An argument must identify which one it explains.
Shaft Recovery Depends on Phase and on the Moving Grip
For a frozen-coefficient linear bending model, write \[ D\ddot\eta+D_d\dot\eta+K_\eta\eta=-C^{\mathsf T}\ddot q. \] With constant symmetric \(D,K_\eta\succ0\) and \(D_d\succeq0\), the modal energy obeys \[ E_\eta=\tfrac12\dot\eta^{\mathsf T}D\dot\eta+\tfrac12\eta^{\mathsf T}K_\eta\eta, \qquad \dot E_\eta=-\dot\eta^{\mathsf T}C^{\mathsf T}\ddot q-\dot\eta^{\mathsf T}D_d\dot\eta. \] The sign of the base acceleration alone does not determine the direction of energy exchange; modal velocity and coupling orientation also enter. In a rotating, deforming model, coefficient derivatives and additional inertial terms must be included before using an energy balance. The frozen equation cannot silently replace the full flexible equations during a large swing.
For one mode driven harmonically at frequency \(\Omega\), the displacement response to force has transfer function \((k-m\Omega^2+i c\Omega)^{-1}\). Stiffness changes amplitude and phase; damping and excitation timing matter as well. The velocity contribution at impact depends on the response at that particular instant. A greater stored elastic energy earlier in the swing does not guarantee a greater final clubhead speed or a preferred face orientation.
MacKenzie and Boucher’s 2016 experiment tested two shafts with similar inertial properties in 33 male golfers, using 14 drives per shaft. The more flexible shaft produced a larger deflection-related speed contribution, accompanied by lower grip angular velocity at impact; the group difference in total clubhead speed was not statistically significant (\(p=0.14\)). Shaft bending and grip orientation also contributed differently to delivered club orientation. The authors considered both mechanical reactions and changed motor patterns as possible explanations. This study supports analyzing the coupled golfer–club system. It does not establish equivalence between shafts, identify one neural mechanism, or guarantee an individual fitting response. Individual significance tests from one session also require attention to multiple comparisons and repeatability. See MacKenzie and Boucher, The Influence of Golf Shaft Stiffness on Grip and Clubhead Kinematics.
The important modeling question is what is held fixed. Prescribing identical grip motion for two shafts estimates an equipment response under a kinematic constraint. Prescribing identical applied forces allows the grip motion to change. Allowing the golfer’s controller to adapt defines a third experiment. Differences among those results are scientifically informative; they should not be collapsed into a single claim that a shaft either adds or does not add speed.
Carry Uncertainty Through the Impact Event
Let \(\Phi(t,t_0)\) propagate infinitesimal state perturbations along a specified smooth trajectory. At a fixed time \(T\), an output \(y=h(x(T),T)\) has initial-state sensitivity \(h_x\Phi(T,t_0)\). Ball contact, however, generally occurs at a state-dependent time defined by \(g(x,t)=0\). Put \(n=\nabla_x g\), let \(F^-\) be the pre-impact vector field, and evaluate quantities at the nominal event. If \(d=g_t+n^{\mathsf T}F^-\ne0\), then \[ \delta t_e=-\frac{n^{\mathsf T}\Phi\delta x_0}{d},\qquad \delta x_e=\left(I-\frac{F^-n^{\mathsf T}}{d}\right)\Phi\delta x_0. \] For an explicitly time-dependent impact output, \(\delta y=h_x\delta x_e+h_t\delta t_e\). These are sensitivities of the state and output evaluated at each perturbed trajectory’s own event. Comparing post-impact trajectories at a common later time instead requires the reset derivative and saltation correction. Near grazing contact, \(d\) can be small and the linear event sensitivity can become large or cease to apply.
An error direction invisible in today’s club-position Jacobian can affect tomorrow’s impact through \(\Phi\). Conversely, reducing every joint error equally can waste effort on directions that have little influence on the chosen outcome. To compare sensitivities, specify scales: if \(\delta x=D_x\delta\bar x\) and \(\delta\bar y=D_y^{-1}\delta y\), use the singular values of \(D_y^{-1}J D_x\), not an unscaled matrix mixing meters, radians, and rates. A singular value is a local worst-case directional gain under those scales, not the expected error of a population of swings.
Parameter uncertainty enters alongside state uncertainty. For \(\delta y=J_x\delta x+J_p\delta p\), propagate both covariance blocks and their cross covariance: \[ \Sigma_y\simeq J_x\Sigma_xJ_x^{\mathsf T}+J_p\Sigma_pJ_p^{\mathsf T} +J_x\Sigma_{xp}J_p^{\mathsf T}+J_p\Sigma_{xp}^{\mathsf T}J_x^{\mathsf T}. \] Mass properties, shaft stiffness, contact location, face orientation, ball properties, and collision assumptions can alter the final prediction. This first-order expression requires small uncertainties and a differentiable local model; it does not cover unmodeled contact changes or large distribution tails. A high-quality golf model should report which uncertainty dominates the output it seeks to explain.
Robotics: Reachability, Redundancy, and a Verifiable Controller
Check the Reference Before Designing Its Controller
For a planar three-link robot with \(l_i=0.30\) m, the end position is \[ p(q)=\sum_{i=1}^3 l_i\begin{bmatrix}\cos\theta_i\\\sin\theta_i\end{bmatrix},\qquad \theta_i=\sum_{j=1}^i q_j. \] A circle centered at \((0.60,0.60)\) m with radius \(0.20\) m is not entirely reachable: its farthest point is about \(1.049\) m from the base, beyond the \(0.90\) m total link length. No feedback gain makes that reference feasible.
Instead choose the manufactured two-second reference \[ p_d(t)=\begin{bmatrix}0.45+0.05\cos(\pi t)\\0.15+0.05\sin(\pi t)\end{bmatrix}\ \mathrm{m}, \qquad \psi_d=\pi/4,\qquad 0\le t\le2\ \mathrm{s}. \] At a quarter period, \(t=0.50\) s, the position is \((0.45,0.20)\) m. Subtract the last link to obtain the wrist target \(w=p_d-0.30(\cos\psi_d,\sin\psi_d)^{\mathsf T}\). Its distance from the base stays between approximately \(0.196\) and \(0.296\) m, strictly within the two-link workspace and away from its folded and extended singularities.
A continuous elbow branch is \[ \begin{aligned} q_{2,d}&=\arccos\left(\frac{w_x^2+w_y^2-0.18}{0.18}\right),\\ q_{1,d}&=\operatorname{atan2}(w_y,w_x) -\operatorname{atan2}(0.30\sin q_{2,d},0.30+0.30\cos q_{2,d}),\\ q_{3,d}&=\psi_d-q_{1,d}-q_{2,d}. \end{aligned} \] Angles must be unwrapped consistently when differentiating the reference. Substituting these angles into forward kinematics verifies the pose; differentiating the resulting smooth branch supplies reference velocities and accelerations. Workspace feasibility alone does not establish torque feasibility, joint-limit clearance, or collision avoidance.
A Position Pullback Has a Nullspace
The position Jacobian has columns \[ J_{p,:,j}=\sum_{i=j}^3 l_i\begin{bmatrix}-\sin\theta_i\\\cos\theta_i\end{bmatrix}. \] It is a \(2\times3\) matrix. At a regular configuration it has rank two, so the position-error pullback \(Q_q=J_p^{\mathsf T}W_pJ_p\), for \(W_p\succ0\), has rank two and a one-dimensional nullspace. It cannot have three positive eigenvalues merely because the robot is away from a position singularity. A determinant of \(J_p\) itself is not defined; manipulability can instead use \(\det(J_pJ_p^{\mathsf T})\) with declared units and scales.
Adding orientation yields the square pose Jacobian with third row \((1,1,1)\). It is nonsingular when the first two links are not collinear. These are different tasks with different nullspaces. In a humanoid swing, the same distinction appears when a cost constrains club position while leaving face orientation or redundant body motion weakly penalized.
Adding \(\epsilon I\) to \(Q_q\) makes that particular coordinate quadratic form positive definite for \(\epsilon>0\). It changes the cost in every direction, not only the nullspace. It does not create a physical barrier against singularities, guarantee a feasible trajectory, or prove contraction. A state metric for a torque-driven three-joint robot also needs six-state treatment; a \(3\times3\) posture penalty alone does not measure velocity error.
An Exact Local Example With Its Assumptions Exposed
For an ideal fully actuated robot with known dynamics \(M(q)\ddot q+h(q,\dot q)=u\), computed torque can impose a desired acceleration: \[ u=M(q)(\ddot q_d+v)+h(q,\dot q),\qquad \ddot e=v,\qquad e=q-q_d. \] Use normalized angle and time coordinates to consider one error channel \(z=(e,e')\), where the prime denotes differentiation with respect to dimensionless time \(\tau=t/t_s\). The physical virtual acceleration is \((q_s/t_s^2)v\) when \(e=(q-q_d)/q_s\). In these normalized variables, \[ z'=A_0z+B_0v,\qquad A_0=\begin{bmatrix}0&1\\0&0\end{bmatrix},\qquad B_0=\begin{bmatrix}0\\1\end{bmatrix}. \] For the infinite-horizon objective \(\int_0^\infty(z^{\mathsf T}z+v^2)\,d\tau\), direct substitution into the continuous Riccati equation gives \[ S=\begin{bmatrix}\sqrt3&1\\1&\sqrt3\end{bmatrix},\qquad K=B_0^{\mathsf T}S=\begin{bmatrix}1&\sqrt3\end{bmatrix},\qquad v=-Kz. \] This standard double-integrator example is also derived in MIT’s Linear Quadratic Regulators notes. The cost is on normalized error and virtual acceleration; it is not a claim about actual motor electrical energy or a human’s physiological objective.
Let \(A_c=A_0-B_0K\) and \(W=I+K^{\mathsf T}K\). The Riccati identity is \[ A_c^{\mathsf T}S+SA_c=-W. \] It follows that the metric length of an error difference contracts at any rate \(\lambda\le\tfrac12\lambda_{\min}(W,S)\) in normalized time, where the generalized eigenvalue uses \(Wv=\mu Sv\). Here the smallest generalized eigenvalue is \(\sqrt3-1/\sqrt2\), so \(\lambda=(\sqrt3-1/\sqrt2)/2\). Physical-time rate is \(\lambda/t_s\). The other common quantity, the real part of a closed-loop eigenvalue, is not interchangeable with this particular metric’s instantaneous bound.
For three independent ideal error channels, repeat this block in a six-state metric after fixing the state ordering. The exact error equation makes the certificate valid in the chosen unwrapped coordinate chart as long as computed torque remains realizable and the model is exact. Saturation, contact changes, uncertain inertia, actuator dynamics, and branch changes invalidate that exact cancellation and require a new residual analysis. The certificate should not be advertised as a verified global property of an arbitrary robot or human.
Finite-Horizon and Sampled Certificates Need Their Own Equations
For a time-varying linear error system with the finite-horizon Riccati solution \(S(t)\) and \(K(t)=R^{-1}B^{\mathsf T}S\), the closed-loop identity is \[ \dot S+A_c^{\mathsf T}S+SA_c=-Q-K^{\mathsf T}RK. \] Uniform positive definiteness of \(S\) and a uniform generalized lower bound on the right-hand quadratic form are needed for a uniform contraction statement. A zero terminal cost can make \(S\) degenerate at the horizon. Plotting a decreasing Riccati trace is neither that test nor a substitute for it.
For sampled error dynamics \(z_{k+1}=A_{c,k}z_k\), the corresponding requirement is \(A_{c,k}^{\mathsf T}S_{k+1}A_{c,k}\preceq\rho^2 S_k\) with \(\rho<1\). In the fixed-metric case, compute \(\rho^2\) as the largest generalized eigenvalue of \((A_c^{\mathsf T}SA_c,S)\). It is not generally the squared spectral radius of \(A_c\). These distinctions matter whenever a controller designed in continuous time is executed on sampled measurements.
Spacecraft: Respect the Orbit and the Meaning of the Cost
Relative Motion About a Circular Chief Orbit
Let \(x\) be radial displacement, \(y\) along-track displacement, and \(z\) normal displacement in the rotating frame of a circular chief orbit. For small relative separation, the Clohessy–Wiltshire linearization is \[ \ddot x=3n^2x+2n\dot y+u_x,\qquad \ddot y=-2n\dot x+u_y,\qquad \ddot z=-n^2z+u_z. \] The inputs are accelerations. The circular-orbit and small-separation assumptions are part of the model; a different axis convention changes the displayed signs. See MIT Astrodynamics, Lecture 26 for the relative-motion derivation.
For radius \(r_0=6.7\times10^6\) m and \(\mu=3.98600435507\times10^{14}\) m\(^3\)/s\(^2\), \(n=\sqrt{\mu/r_0^3}\simeq0.001151\) rad/s and \(T=2\pi/n\simeq5458\) s, about \(90.96\) minutes. This is a low Earth orbit, not a geosynchronous orbit. The Earth gravitational parameter is from JPL’s DE440 astrodynamic parameters.
For \(X=(x,y,z,\dot x,\dot y,\dot z)\), write \(\dot X=A X+B u\) with \[ A=\begin{bmatrix}0&I_3\\N&C_o\end{bmatrix},\quad B=\begin{bmatrix}0\\I_3\end{bmatrix},\quad N=\operatorname{diag}(3n^2,0,-n^2),\quad C_o=\begin{bmatrix}0&2n&0\\-2n&0&0\\0&0&0\end{bmatrix}. \] The blocks act on position and velocity with different units; it is misleading to attach a single dimensional unit to every entry of this unscaled state matrix.
A reproducible reference can be manufactured rather than claimed as flight validation. Over \(t_f=900\) s, set \(s=t/t_f\), \(p(s)=10s^3-15s^4+6s^5\), and \(y_d=1000-990p(s)\) m, with \(x_d=z_d=0\). Then \(u_{x,d}=-2n\dot y_d\), \(u_{y,d}=\ddot y_d\), and \(u_{z,d}=0\) satisfy the linear equations. The reference starts and ends at rest, with zero endpoint acceleration, and ends at a 10 m along-track standoff. It is not a completed docking maneuver. Thruster limits, approach corridors, collision geometry, navigation error, and nonlinear orbital-model error still require evaluation.
Quadratic Effort Is Not Propellant Consumption
An LQR objective \(\int(e^{\mathsf T}Qe+\delta u^{\mathsf T}R\delta u)\,dt\) balances weighted tracking error and squared control effort around a feasible reference. Position and velocity weights require scales; \(Q\) here is a cost weight, not a process-noise covariance. A terminal penalty is a soft cost, not an exact terminal docking constraint.
For an ideal thruster with constant specific impulse, \(\dot m=-\|F\|/(I_{sp}g_0)\) describes propellant use. The actual thruster allocation and operating model can require modifications. With \(F=mu\), mass changes can enter the dynamics. The time integral of squared acceleration is therefore not a direct fuel measurement. A fuel-optimal or fuel-limited claim requires the appropriate objective, thrust model, mass evolution, and constraints. Neither a quadratic regulator nor an unreported simulation establishes a numerical percentage of fuel savings.
To implement a sampled controller under zero-order-held acceleration, compute \[ \exp\left(\begin{bmatrix}A&B\\0&0\end{bmatrix}\Delta t\right) =\begin{bmatrix}A_d&B_d\\0&I_3\end{bmatrix}. \] Design a discrete regulator for \((A_d,B_d)\) and evaluate \(A_d-B_dK_d\). Subtracting a discrete gain from the continuous matrix and interpreting the result as a sampled stability certificate mixes two different systems. A covariance propagation or Monte Carlo study can then test navigation uncertainty, but it must report its noise model and cannot validate unmodeled thruster or contact physics.
Vehicles: Input Coordinates and Speed Change the Argument
An Exact Affine Bicycle Model
For the kinematic bicycle with rear-axle position \((x,y)\), heading \(\psi\), speed \(v\), wheelbase \(L\), steering angle \(\delta\), and longitudinal acceleration \(a\), \[ \dot x=v\cos\psi,\quad \dot y=v\sin\psi,\quad \dot\psi=\frac{v}{L}\tan\delta,\quad \dot v=a. \] This model is nonlinear in the steering-angle input. With curvature \(\kappa=\tan\delta/L\), it is exactly control affine: \[ \dot X=\begin{bmatrix}v\cos\psi\\v\sin\psi\\0\\0\end{bmatrix} +\begin{bmatrix}0\\0\\0\\1\end{bmatrix}a +\begin{bmatrix}0\\0\\v\\0\end{bmatrix}\kappa. \] Alternatively, append steering angle to the state and use steering rate as an input. A Taylor expansion in \(\delta\) around a reference gives a local input approximation, not the exact original affine model. Steering limits transform to curvature limits, and steering-rate constraints do not disappear under this reparameterization.
The no-side-slip condition restricts velocities. It does not by itself reduce the dimension of the configuration space: sequences of admissible motions can produce displacements not available instantaneously. Lie brackets describe this local mechanism, as illustrated in Modern Robotics, Controllability of Wheeled Mobile Robots. A rank calculation alone should not be turned into an unrestricted global parking guarantee; reversibility, admissible controls, obstacles, and the particular controllability theorem’s assumptions matter.
A Speed Schedule Requires Its Own Dynamics
For a straight reference and small heading error, with fixed \(v>0\), the lateral error equations are \(\dot e_y=v e_\psi\) and \(\dot e_\psi=v\kappa\). The curvature feedback \[ \kappa=-\frac{\omega_n^2}{v^2}e_y-\frac{2\zeta\omega_n}{v}e_\psi \] gives \(\ddot e_y+2\zeta\omega_n\dot e_y+\omega_n^2e_y=0\). For positive \(\omega_n\) and \(\zeta\), this linearized constant-speed model is stable. The position and heading gains decrease as speed increases for this fixed target polynomial. There is no general rule that all gains should increase with speed.
If \(v\) varies, differentiation introduces \((\dot v/v)\dot e_y\), and simply inserting the current speed into the same law no longer gives the stated constant-coefficient equation. The gain formula is singular at zero speed. Reverse motion, curved references, steering dynamics, and saturation also need separate treatment. Frozen-speed eigenvalues do not certify stability under arbitrary switching or rapid scheduling.
At fixed steering angle, yaw rate scales as \(v\) and the no-slip lateral acceleration as \(v^2\tan\delta/L\). For \(L=3\) m, \(\delta=0.10\) rad, and \(v=25\) m/s, the latter is about \(20.9\) m/s\(^2\). That requested acceleration signals that the simple kinematic assumptions need scrutiny; the equation does not establish that the tires can supply it. A dynamic tire model and friction constraints determine feasibility. Comparing unscaled norms that mix position rates, heading rate, and longitudinal acceleration cannot establish control dominance at high speed.
The connection to golf is methodological. A parameter change can alter the dynamics, input map, feasible controls, and error metric simultaneously. A gain schedule or equipment change must be evaluated within that coupled system; conclusions from a fixed-parameter calculation do not automatically survive the change.
Reproducible Computation and Evidence
Check an Equation Before Interpreting Its Output
The following complete Python example checks the normalized double-integrator result. It requires NumPy and SciPy. The residual verifies the equation solved; the generalized eigenvalue verifies the stated metric inequality. Neither check validates an unmodeled actuator.
import numpy as np
from scipy.linalg import eigvalsh, solve_continuous_are
A = np.array([[0.0, 1.0], [0.0, 0.0]])
B = np.array([[0.0], [1.0]])
Q = np.eye(2)
R = np.eye(1)
S = solve_continuous_are(A, B, Q, R)
K = np.linalg.solve(R, B.T @ S)
Ac = A - B @ K
W = Q + K.T @ R @ K
residual = Ac.T @ S + S @ Ac + W
assert np.linalg.norm(residual) < 1e-10
rate = 0.5 * eigvalsh(W, S)[0]
assert np.isclose(rate, (np.sqrt(3.0) - 1.0 / np.sqrt(2.0)) / 2.0)For an equilibrium regulator, the physical command includes any nonzero equilibrium input. For a trajectory regulator, it includes a dynamically feasible reference input and the tracking correction. Omitting that feedforward term changes the error dynamics. Apply limits and actuator dynamics before interpreting simulated performance; clipping an ideal command invalidates an unconstrained certificate unless the resulting residual is bounded in a new analysis.
In MATLAB, dlqr(A,B,Q,R,N) solves an infinite-horizon discrete LQR problem; the fifth argument is the state–input cross-weight matrix \(N\), not a horizon length. A finite-horizon problem requires a backward recursion with a declared terminal cost and time-indexed gains. The MathWorks dlqr documentation specifies the cross term and discrete Riccati equation. Likewise, an outline containing undefined model objects or placeholder matrices is an implementation sketch, not an executable validation.
Separate Verification, Calibration, and Validation
Verification asks whether the equations were implemented and solved correctly. For the golf calculation this includes Cartesian-versus-joint kinetic energy, mass symmetry and definiteness under admissible parameters, correct signed mode normalization, correct Schur elimination, unforced energy conservation, and step refinement. A reproducible table is valuable only after these model identities are satisfied.
Calibration estimates model parameters from data. Validation compares held-out predictions with observations for a declared use. Fitting the same swing used to assess accuracy is not an independent prediction test. A claimed error percentage needs the measured quantity, normalization, sample, uncertainty, and protocol; it cannot be inferred from the apparent plausibility of a simulated trajectory.
For perturbation experiments, align the model and measurement definitions. Separate fixed-time errors from event-time errors, and choose whether the comparison is between two closed-loop trajectories or between a trajectory and its nominal reference. Measurement noise can dominate a ratio with a small initial-error denominator. A worst-case contraction bound need not be attained by typical perturbations, and a confidence interval containing a predicted bound is not a validation criterion by itself. Test the inequality with appropriate uncertainty, check the assumed region and disturbances, and report violations as well as agreement.
For the golf swing, an informative progression is to verify rigid-body energy accounting, add and identify shaft dynamics, verify contact and impact behavior, then assess the coupled model on held-out swings and equipment conditions. Measure grip and head motion together when studying shaft effects. Include force or moment information when testing an energy-transfer explanation. Include event timing and output uncertainty when testing an accuracy explanation. Each additional model component should resolve a stated discrepancy or prediction need rather than merely increase the number of states.
Exercises That Test the Argument
- Derive the six independent entries of \(M_r\) from Cartesian center-of-mass velocities. Check at random postures and velocities that both kinetic-energy calculations agree. Explain why an endpoint singularity need not make \(M_r\) singular.
- Assemble the beam’s three inertia blocks from the same material-point velocity. Repeat with twice the beam mass. Show which rigid-coordinate diagonal term gains \(m_s l_3^2/3\), and test the complete kinetic energy rather than only a matrix determinant.
- Change the sign of the second mode and transform its amplitude and coupling consistently. Show that physical displacement and kinetic energy are unchanged. Explain why assigning two arbitrary frequencies does not identify a uniform beam’s \(EI\).
- Verify the two block acceleration equations by substituting the Schur solution. Include a nonzero modal bias so an omitted elastic term cannot pass the check. Distinguish this test from integrating the rigid six-state counterfactual.
- Refine the rigid zero-torque integration step and report both endpoint convergence and mechanical-energy error. Explain why the displayed acceleration norm is not a finite drift-control ratio or an estimate of muscular effort.
- Construct a motion with decreasing proximal angular speed but no positive distal energy transfer across the chosen joint. State the external powers and show the segment energy balance. Identify the extra measurements required to interpret a human speed-peak sequence.
- For two shaft stiffnesses, compare identical prescribed grip motion, identical applied force histories, and an adapting controller. State which causal question each comparison addresses and which output would distinguish the hypotheses.
- Reconstruct the robot circle using the given inverse kinematics. Compute the nullspace of \(J_p\) and compare position and pose pullbacks. Test whether adding \(\epsilon I\) imposes any obstacle or joint-limit constraint.
- Verify the exact Riccati solution and its generalized metric rate. Add a torque limit to the computed-torque model and identify where the ideal proof ceases to apply rather than extending its rate by assertion.
- Substitute the quintic spacecraft reference into the relative-motion equations. Compare integrated squared acceleration with a declared propellant model. Report the assumptions needed before calling the final standoff a docking solution.
- Derive the additional lateral-error term for variable vehicle speed. Test a rapid speed change and a steering limit. Explain why frozen-speed eigenvalues and kinematic lateral acceleration do not establish dynamic feasibility.
- For a transverse impact event, compare finite-difference sensitivities with the event-time formula and a fixed-time formula. Repeat near grazing. Propagate parameter as well as state uncertainty, using declared input and output scales.
Citation
@online{untitled,
author = {},
title = {Applications: {Biomechanics,} {Robotics,} and {Beyond}},
url = {https://affinedrift.com/articles/The_Geometry_of_Motion/quarto/ch08_applications.html},
langid = {en}
}