Contraction Theory Meets Tangent Spaces: Stability, Optimality, and Scope Limits

What a Differential Stability Certificate Establishes

Derivations and counterexamples connecting differential stability, Riccati value functions, constrained control, and measurable golf-swing outcomes.

Two copies of a motion can start slightly apart. Contraction theory asks whether that gap shrinks under a given model and feedback rule. This article links that test to optimal control and to golf. For a longer plain-language walk-through, read the accessible explanation.

One Calculation, Two Jobs

Under stated conditions, the same matrix that prices a mistake can also prove that nearby mistakes shrink. Neither optimality nor a positive matrix alone gives that proof.

Limits of the Proof

The result holds for one model, one feedback rule, one way of measuring error, and one region. A rubber band is a loose picture: a spring alone does not remove the gap, and the metric is not a real force.

What Golf Needs

A repeatable swing need not erase every difference. Some variation may not change club delivery, while a small error elsewhere may. Those golf questions need measurement, and the math does not settle them.

Key Takeaway: A feedback stability certificate describes a declared model and region. It does not show how any golfer swings.

Abstract

Contraction analysis asks whether neighboring solutions approach one another under a specified dynamics model and input policy. Trajectory optimization asks which admissible motion minimizes a declared objective. Their common use of differential dynamics creates a useful connection: under explicit definiteness and boundedness conditions, a linear-quadratic value matrix also certifies contraction of the associated closed loop. Neither optimality nor a positive matrix alone supplies that certificate.

This article derives the connection, distinguishes continuous and discrete time, and develops counterexamples to stronger interpretations. It then treats metric geometry, numerical verification, musculoskeletal actuation, stochastic variability, task-space control, and contact events. The golf application connects geometry and stored energy to feasible corrections and the delivered impact state. Its empirical questions remain separate from the mathematical identities.

1 Reading the Physical Argument

A repeatable swing need not erase every state difference. It may tolerate variation in directions that have little influence on club delivery while controlling a small number of consequential errors. A contraction metric is one way to test a stronger model claim: all differential errors measured by that metric decrease throughout a specified region and time interval.

The metric changes how error is measured; feedback and mechanics change how error evolves. Calling a metric curved does not make it a physical restoring force. Earlier work, geometry, inertia, grip compliance, activation history, and ground contact determine the state and the available responses. A stability calculation must include those mechanisms in its dynamics before drawing conclusions about their effects.

2 Differential Dynamics and the Correct Metric Test

Let the smooth closed-loop dynamics be \[ \dot x=F(x,t)=f(x,\pi(x,t),t),\qquad A_{\mathrm{cl}}=D_xF. \tag{1}\] Here \(\pi\) is the actual input policy. Differentiating while holding time fixed gives the exact variational equation \[ \delta\dot x=A_{\mathrm{cl}}(x,t)\delta x. \tag{2}\] For an open-loop plant with independently perturbed input, the equation is instead \(\delta\dot x=f_x\delta x+f_u\delta u\). Closing the loop adds \(\delta u=\pi_x\delta x\), so \(A_{\mathrm{cl}}=f_x+f_u\pi_x\). Differentiating the plant with input held fixed does not test feedback stability.

Choose a smooth symmetric metric \(M(x,t)\) and define differential squared length \(V=\delta x^\top M\delta x\). Its derivative is \[ \begin{aligned} \dot V&=\delta x^\top\mathcal H\delta x,\\ \mathcal H&=A_{\mathrm{cl}}^\top M+MA_{\mathrm{cl}}+\dot M,\\ \dot M&=\partial_tM+\sum_iF_i\partial_{x_i}M. \end{aligned} \tag{3}\] The material derivative uses the same closed-loop flow. A sequence of matrices evaluated along one trajectory does not define the spatial derivatives of a metric in its neighborhood.

The contraction inequality with length-decay rate \(\lambda>0\) is \[ \mathcal H\preceq-2\lambda M. \tag{4}\] Its normalized symmetric matrix is \[ M^{-1/2}\mathcal H M^{-1/2}\preceq-2\lambda I. \tag{5}\] The multiplication order matters. With \[ A=\begin{bmatrix}-1&3\\0&-2\end{bmatrix},\qquad M=\begin{bmatrix}1/2&1/2\\1/2&1\end{bmatrix}, \] the correct \(A^\top M+MA=-I\). The different expression \(AM+MA^\top\) has a positive eigenvalue approximately \(2.354\) and tests a different object.

For a constant LTI matrix and constant metric this specializes to \[ A^\top M+MA\preceq-2\lambda M,\qquad M\succ0. \tag{6}\]

If \(M=\Theta^\top\Theta\) with invertible \(\Theta\), differential components \(\delta z=\Theta\delta x\) obey \[ \delta\dot z=\mathcal F\delta z,\qquad \mathcal F=(\dot\Theta+\Theta A_{\mathrm{cl}})\Theta^{-1}. \tag{7}\] Then \(\operatorname{sym}\mathcal F\preceq-\lambda I\) is equivalent to the metric inequality. The differential transformation need not integrate to a global coordinate map \(z(x)\).

2.1 From Differential Length to Finite Separation

Integrating the differential inequality gives \[ \|\delta x(t)\|_{M(x(t),t)} \leq e^{-\lambda(t-t_0)}\|\delta x(t_0)\|_{M(x(t_0),t_0)}. \tag{8}\] For finite separation, connect the initial states by a smooth path and evolve every point of that path under the same flow. If this path of solutions remains inside the certified region, its metric length shrinks. Taking the infimum over admissible connecting paths yields a corresponding bound on Riemannian distance. A finite difference \(x_2-x_1\) cannot simply be substituted for \(\delta x\) in the exact differential equation of a nonlinear system.

State existence, path containment, and metric bounds belong to the conclusion. Uniform bounds \(0<m_-I\preceq M\preceq m_+I\) permit comparison with Euclidean errors, including the conditioning factor \(\sqrt{m_+/m_-}\). An invariant region or an invariant tube about a reference can provide the required containment. A condition checked only at the reference does not establish the tube. These distinctions follow the path-based formulation in Lohmiller and Slotine.

Contraction shrinks every differential length in the chosen metric. Volume shrinkage is weaker: \(\dot x=\operatorname{diag}(1,-2)x\) shrinks area while one direction expands. Convergence between solutions also does not imply bounded motion: \(\dot x=-x+e^{2t}\) has contracting differences while its solutions grow without bound. For an autonomous system, strict full-state contraction cannot support a nonconstant periodic orbit in the same certified domain, because phase-shifted solutions would have to converge at equal time. Orbital or transverse contraction is a different statement.

2.2 Disturbances and Interconnections

Suppose one trajectory receives an added disturbance whose norm in the certified metric is bounded by \(\bar d\). Under the same region and regularity conditions, a comparison inequality gives \[ d_M(t)\leq e^{-\lambda(t-t_0)}d_M(t_0) +\frac{\bar d}{\lambda}\left(1-e^{-\lambda(t-t_0)}\right). \] Persistent disturbance generally produces a nonzero error bound. It does not guarantee convergence to the undisturbed trajectory, and the perturbed trajectory must remain in the region where the bound applies.

Independent contracting subsystems driven by the same prescribed signal admit a block-diagonal contraction metric. Mechanically coupled joints are not independent subsystems. Feedback interconnection adds off-diagonal Jacobian blocks whose symmetric contribution must be bounded. A bounded hierarchical coupling can sometimes be handled by rescaling a composite metric; arbitrary mutual coupling cannot be ignored. For example, two isolated equations \(\dot x_i=-x_i\) are stable, but coupling them as \(\dot x_1=-x_1+2x_2\), \(\dot x_2=2x_1-x_2\) creates an expanding eigenvalue \(1\).

3 Riccati Value Functions and Contraction

3.1 The Continuous-Time Identity

Consider a linear time-varying perturbation system \(\dot\xi=A(t)\xi+B(t)v\) and the quadratic cost \[ J=\xi(T)^\top S_T\xi(T) +\int_{t_0}^T\left(\xi^\top Q\xi+2\xi^\top Nv+v^\top Rv\right)dt, \qquad R\succ0. \] Assume the cost block is positive semidefinite and the Riccati solution exists. Completing the square in dynamic programming gives \[ -\dot S=A^\top S+SA-(SB+N)R^{-1}(B^\top S+N^\top)+Q, \qquad K=R^{-1}(B^\top S+N^\top). \tag{9}\] For the minimizing policy \(v=-K\xi\), let \(A_c=A-BK\) and \[ W=Q-NK-K^\top N^\top+K^\top RK. \] Substitution yields the exact identity \[ \dot S+A_c^\top S+SA_c=-W. \tag{10}\] Thus \(S\) is a contraction metric on an interval if it is uniformly positive definite and bounded there and \(W\succeq2\lambda S\) holds uniformly. For zero cross cost, \(W=Q+K^\top RK\). The open-loop ARE instead gives \(A^\top S+SA=SBR^{-1}B^\top S-Q\); the negative feedback contribution appears only after forming \(A_c\).

For infinite-horizon LTI LQR with \(N=0\), standard stabilizability and detectability assumptions with \(Q\succeq0\), \(R\succ0\) yield a stabilizing positive-semidefinite solution. Positive definiteness and a strict decay bound in that particular matrix require their own checks. For example, \(A=-1\), \(B=1\), \(Q=0\), \(R=1\) has optimal cost \(S=0\) and a stable uncontrolled plant, but zero is not a Riemannian metric. With \(Q\succ0\) the usual stabilizing solution is positive definite. The MIT LQR treatment develops the finite- and infinite-horizon optimization problems.

3.2 Finite-Horizon Decrease Is Not an Asymptotic Certificate

A terminal weight may vanish. In the scalar uncontrolled system \(\dot\xi=\xi\) with running cost \(\xi^2\) and zero terminal cost, \[ S(t)=\frac{e^{2(T-t)}-1}{2},\qquad \dot S+2S=-1. \] The value \(S(t)\xi(t)^2\) decreases as the remaining horizon disappears, while \(|\xi(t)|\) grows. At \(T\), \(S(T)=0\). There is no uniformly positive metric through the endpoint and no claim about behavior after the horizon. This example isolates why value decrease and physical error decay must be distinguished.

In a nonlinear tracking problem, linearization along a nominal trajectory supplies a variational candidate. To extend its result to a finite neighborhood, verify the closed-loop derivative and metric inequality throughout a tube, with strict margin, model residual bounds, feasible controls, and containment. Recomputing LQR at sample points does not supply these conditions automatically.

3.3 Discrete-Time Dynamics Require a Different Test

For a discrete closed loop \(x_{k+1}=F_k(x_k)\) with differential map \(A_{c,k}=D F_k\), the contraction condition is \[ A_{c,k}^\top M_{k+1}(F_k(x))A_{c,k}\preceq\rho M_k(x),\qquad0<\rho<1. \tag{11}\] The metric is evaluated at the next state on the left. At sample interval \(\Delta t\), the length-decay rate is \(-\log(\rho)/(2\Delta t)\). A continuous Lyapunov inequality applied to the Jacobian of a discrete integrator is not this condition.

For zero-cross-cost discrete LQR, the Riccati identity is \[ (A_k-B_kK_k)^\top S_{k+1}(A_k-B_kK_k) =S_k-Q_k-K_k^\top R_kK_k. \tag{12}\] If \(S_k\) has uniform positive bounds and \(Q_k+K_k^\top R_kK_k\succeq(1-\rho)S_k\), it supplies the desired contraction inequality. The normalized matrix inequality, not the spectral radius of each update, establishes the one-step rate.

For example, \[ F_1=\begin{bmatrix}0.5&2\\0&0.5\end{bmatrix},\qquad F_2=F_1^\top. \] Both matrices have eigenvalues \(0.5\), yet \(F_2F_1\) has an eigenvalue approximately \(4.486\). Alternation amplifies errors. Time-varying stability depends on products or a valid changing metric, not independent eigenvalue tests. Even in an LTI system, nonnormal transient amplification separates an asymptotic spectral rate from monotone decay in a chosen norm.

3.4 Worked Double-Integrator Certificate

Take \[ A=\begin{bmatrix}0&1\\0&0\end{bmatrix},\quad B=\begin{bmatrix}0\\1\end{bmatrix},\quad Q=\operatorname{diag}(16,1),\quad R=1. \] The stabilizing ARE solution and feedback are \[ S=\begin{bmatrix}12&4\\4&3\end{bmatrix},\qquad K=\begin{bmatrix}4&3\end{bmatrix}. \tag{13}\] They give eigenvalues \(-1.5\pm(\sqrt7/2)i\). In the specific metric \(S\), however, the largest certified uniform length-decay rate is \[ \lambda_S=\frac12\lambda_{\min}\left(S^{-1/2}(Q+K^\top K)S^{-1/2}\right) \approx1.276393\ \mathrm{s}^{-1}. \] The other normalized half-eigenvalue is approximately \(1.723607\). Neither the eigenvalue real part \(1.5\) nor a requested design value \(2\) can replace this calculation. Rates of exponential decay have units inverse seconds; they are not oscillation frequencies measured in cycles per second.

LQR minimizes its declared cost, not contraction rate over all feedback policies. Changing the cost weights changes that compromise; a scalar weight ratio cannot prescribe arbitrary multivariable poles. To design a rate directly for an unconstrained LTI plant, set \(X=M^{-1}\) and \(Y=KX\). Congruence transforms the state-feedback inequality into \[ AX+XA^\top-BY-Y^\top B^\top+2\lambda X\preceq0, \qquad X\succ0. \] For fixed \(\lambda\), this is linear in \(X,Y\). Add normalization and any required actuator or region constraints, and recover \(K=YX^{-1}\). This feasibility problem is not automatically an inverse-LQR solution, and imposing input bounds over a state region requires additional constraints.

4 What the Geometry Does and Does Not Say

4.1 Metric Tensors and Coordinates

For a diffeomorphism \(z=\phi(x)\), write \(T=D\phi\). Since \(\delta z=T\delta x\), equal differential lengths require \[ M_z=T^{-\top}M_xT^{-1}. \] The expression \(J^\top M_yJ\) is instead the pullback of a metric from a target through a map \(y=h(x)\) with \(J=Dh\). These are compatible tensor rules in opposite directions. Coordinate invariance of the inequality does not make arbitrary choices of metric equivalent or make their Euclidean conditioning identical.

The bilinear metric and task pullback can be written explicitly as \[ g_x(\xi,\zeta)=\xi^\top M_x\zeta, \tag{14}\] \[ h^*M_y=J^\top M_yJ. \tag{15}\]

A Riccati matrix depending on time alone has no spatial derivatives at each fixed time. Its spatial Riemann curvature is zero in these coordinates. LQR trajectories are solutions of controlled dynamics and an optimization problem; they are not generally geodesics of that matrix. For instance, \(\dot x=-x\) contracts in the flat Euclidean metric, without curvature or a funnel-shaped physical surface.

For a metric \(g\), the Levi–Civita coefficients are \[ \Gamma^i_{jk}=\tfrac12g^{i\ell}(\partial_jg_{\ell k}+\partial_kg_{\ell j}-\partial_\ell g_{jk}). \tag{16}\] They vanish for a spatially constant metric in these coordinates. This fact is compatible with a nonzero contraction rate of a dissipative controlled flow.

The variational equation in coordinates uses ordinary time derivatives. With a chosen torsion-free connection it can be written \(\nabla_t\delta x=\nabla_{\delta x}F\); replacing only the left derivative with a covariant derivative while keeping an ordinary Jacobian on the right omits connection terms. This geometric rewriting does not turn the dynamics into an Euler–Lagrange optimality condition. Differentiating an ODE and varying a dynamics-constrained objective are separate operations.

In components, the covariant identity is \[ \frac{D\delta x^i}{dt}=\left(\partial_jF^i+\Gamma^i_{kj}F^k\right)\delta x^j. \tag{17}\]

4.2 Task Metrics, Inertia, and Stiffness

For a task \(y=h(q)\), \(J^\top W_yJ\) measures local task displacement. It is positive definite only if \(J\) has full column rank when \(W_y\succ0\). Redundancy or singularity can leave task-invisible directions. With a straight two-link unit-length arm, \[ J=\begin{bmatrix}0&0\\2&1\end{bmatrix},\qquad J^\top J=\begin{bmatrix}4&2\\2&1\end{bmatrix}, \] which is singular. It equals the generalized kinetic metric of a unit point mass at the endpoint with massless links. It is not the mass matrix of two unit-mass links. For two uniform unit rods the actual inertia is \(\left[\begin{smallmatrix}8/3&5/6\\5/6&1/3\end{smallmatrix}\right]\), which is positive definite at this configuration.

Mechanical inertia acts on configuration velocities. A contraction metric for state \((q,\dot q)\) acts on both position and velocity perturbations and generally includes cross blocks. Stiffness acts between displacement and force; it is not automatically a full-state contraction metric. Likewise, \(BR^{-1}B^\top\) is a contravariant input map that is often singular in underactuated mechanics, not automatically a positive-definite covariant metric. Fisher information and value Hessians can also be degenerate. Every proposed metric needs a domain, units or nondimensionalization, rank check, and differential inequality.

Minimizing \(\operatorname{tr}M\) without a fixed lower bound is degenerate because the homogeneous contraction inequality permits \(\epsilon M\) for arbitrarily small \(\epsilon>0\). Scaling the ruler does not enlarge a physical basin. With \(I\preceq M\preceq\kappa I\), a fixed linear closed loop and fixed rate produce a meaningful semidefinite feasibility problem. A nonlinear metric search must specify a function class, derivative constraints, a domain, and how a continuum inequality is verified. Solver complexity depends on those choices; no universal \(O(n^4)\) runtime follows from naming the problem an SDP.

For the fixed constant closed loop, the normalized feasibility constraints are \[ I\preceq M\preceq\kappa I,\qquad A_c^\top M+MA_c+2\lambda M\preceq0. \tag{18}\]

5 Combining Trajectory Optimization With Stability Verification

5.1 What DDP Actually Expands

For a discrete dynamics map \(x^+=f_d(x,u)\) and running cost \(\ell\), define the local action value \(\mathcal Q=\ell+V^+(f_d(x,u))\). Along a feasible nominal trajectory, write \(A=f_{d,x}\), \(B=f_{d,u}\), \(p=V_x^+\), and \(P=V_{xx}^+\). The second derivatives include \[ \begin{aligned} \mathcal Q_{xx}&=\ell_{xx}+A^\top PA+\sum_i p_i f_{d,xx}^{,i},\\ \mathcal Q_{ux}&=\ell_{ux}+B^\top PA+\sum_i p_i f_{d,ux}^{,i},\\ \mathcal Q_{uu}&=\ell_{uu}+B^\top PB+\sum_i p_i f_{d,uu}^{,i}. \end{aligned} \] DDP retains these dynamics-Hessian contractions. The iLQR approximation omits them while retaining the relevant cost Hessians and cross terms. In either case, a general optimization step includes a feedforward correction as well as feedback: \[ \delta u=k+L\delta x,\qquad k=-\mathcal Q_{uu}^{-1}\mathcal Q_u,\quad L=-\mathcal Q_{uu}^{-1}\mathcal Q_{ux}, \] when the local control Hessian is positive definite and the step is unconstrained. Regularization, input bounds, line search, and defects in an infeasible nominal rollout change the subproblem. The purely homogeneous LQR recursion describes a special case, not a complete nonlinear optimizer. See MIT trajectory optimization.

A value Hessian can be indefinite away from a suitable local minimum. Even a positive Hessian does not certify the resulting nonlinear feedback policy. Optimization convergence, trajectory feasibility, stability, and empirical model validity are distinct checks.

5.2 A Defensible Contraction Penalty

For a declared metric and continuous closed-loop policy, define the normalized defect \[ d(x,t)=\lambda_{\max}\left(M^{-1/2}\mathcal H M^{-1/2}\right)+2\lambda, \qquad c(x,t)=\max(0,d(x,t))^2. \] Zero penalty means the tested point satisfies the inequality at the requested rate. A finite penalty weight in an objective does not force zero violation. Sampled zero violations also do not certify an unsampled region. The largest-eigenvalue function is nonsmooth at repeated active eigenvalues, so an implementation must specify how derivatives or nonsmooth steps are handled.

If a scalar penalty is added to a cost, its first derivatives modify the cost gradients and its second derivatives modify the cost Hessians, including mixed terms. Adding a gradient vector to a Hessian matrix is dimensionally invalid. If the policy or metric is optimized too, account for their dependence in those derivatives. Use the continuous condition for continuous Jacobians and the next-state metric inequality for discrete maps.

For a candidate policy and metric, a useful verification sequence is:

1. Check the nominal trajectory against the dynamics, contact model, and actuator limits.
2. Define the metric throughout a proposed tube, with explicit positive lower and upper bounds.
3. Differentiate the implemented closed loop, including saturation or its smooth region, activation states, and metric material derivatives.
4. Verify the contraction inequality with numerical margin throughout that tube, and establish trajectory containment. A certified Lipschitz bound can turn a sampled margin into a neighborhood bound; an arbitrary sample spacing cannot.
5. Bound model errors and disturbances, inspect contact-event sensitivities, and report the resulting finite-error or output bound.

For example, if the largest eigenvalue of the normalized defect matrix is at most \(-\eta\) at a point and a proven Lipschitz constant is \(L\), a ball of radius \(r<\eta/L\) retains negativity, provided the metric and policy remain regular there. This verifies a local inequality; an invariant tube and time-uniform bounds still need proof. A finite-horizon capture set under one policy is not automatically an asymptotic basin of attraction or a control-reachable set under all policies.

Control contraction metrics provide a separate constructive framework for differential feedback design. The original Manchester–Slotine paper and its integrability amendment specify technical conditions for turning differential controls into finite controls. A pointwise Riccati matrix or a soft penalty should not be presented as an implementation of those theorems without checking their hypotheses.

6 Musculoskeletal Mechanics and Variability

6.1 Force, Activation, and Full Mechanical State

Moment arms map muscle forces to joint torque: \[ \tau=R(q)F_m(q,\dot q,a,z)+\tau_{\mathrm{other}}. \] Here \(a\) denotes activation and \(z\) any additional tendon, fatigue, or history state. Multiplying a moment-arm matrix, measured in length, by dimensionless activation alone does not produce torque. A simplified constitutive law can be affine in activation at fixed muscle state; its force-length, force-velocity, tendon, and excitation dynamics must still be specified. The OpenSim activation model illustrates why excitation and activation occupy different roles.

A low-dimensional factorization \(a\approx Wc\) is a candidate description of recorded activation patterns. Its columns live in muscle-activation space, not directly in the tangent space of the mechanical state. The resulting input directions involve the muscle model, moment arms, inertia, and constraints. A low-rank factorization does not establish that the nervous system selects \(W\) to maximize contraction, or even that the factors identify unique neural modules.

For mechanics \[ D(q)\ddot q+h(q,v)=B(q)\tau,\qquad v=\dot q, \] define \(a_q=D^{-1}(B\tau-h)\) as the acceleration function. Its configuration derivative in direction \(\delta q\) is \[ D_qa_q[\delta q]=D^{-1}\left(D_qB[\delta q]\tau-D_qh[\delta q]-D_qD[\delta q]a_q\right), \] and \(D_va_q=-D^{-1}D_vh\), \(D_\tau a_q=D^{-1}B\) under the stated \(q\)-only inertia and input map. The first-order state Jacobian has top blocks \((0,I)\) and bottom blocks \((D_qa_q,D_va_q)\). Omitting the derivative of inertia or the velocity dependence of the Coriolis load changes the predicted perturbation dynamics. Closed loops add reaction derivatives; muscle-driven models add their internal-state blocks.

6.2 Impedance Is Part of the Dynamics

A local endpoint law \(\delta F=-K_y\delta y-C_y\delta\dot y\) contributes stiffness and damping. The simple pullback \(J^\top K_yJ\) is the unloaded linearization at a matching reference; differentiating \(J^\top F\) at nonzero preload also produces a geometric term involving \(D_qJ\). Along motion, velocity and Jacobian-rate terms matter as well.

Finite stiffness can redirect an error without eliminating a degree of freedom. A spring without damping stores and returns energy; it does not make trajectories converge. A contraction test therefore belongs to the full position/velocity/internal-state dynamics. Co-contraction may alter impedance and later force production without changing instantaneous net torque, and it can have metabolic cost. Neither an inferred torque nor a metric identifies that cost uniquely.

6.3 Stochastic Errors Need a Noise Model

For fully observed additive-noise LQ control, \[ dx=(Ax+Bu)dt+G\,dW_t,\qquad \mathbb E[dW_t\,dW_t^\top]=I\,dt, \] the expected quadratic-cost problem has the deterministic Riccati equation for its state-quadratic part. Noise contributes an additive value term involving \(\operatorname{tr}(SGG^\top)\). Under \(u=-Kx\), the covariance evolves with the closed-loop matrix: \[ \dot\Sigma_x=A_c\Sigma_x+\Sigma_xA_c^\top+GG^\top. \] For \(dx=-x\,dt+\sigma\,dW_t\) from a known initial state, the variance is \(\sigma^2(1-e^{-2t})/2\). Mean error decays, but variance approaches a nonzero floor. Two trajectories driven by exactly the same additive noise have deterministic difference dynamics; independent noise does not cancel. These are different notions of stochastic contraction and synchronization.

Risk-sensitive exponential costs can change the Riccati equation, but their noise-dependent term follows from the specified criterion, normalization, and admissibility conditions. It should not be inserted into ordinary expected-quadratic LQR by analogy. Additive-noise covariance alone also does not establish that a human must always increase effort to maintain a particular outcome; the objective, feedback, signal-dependent noise, and feasible muscle responses need a model and data.

7 Task-Space Tracking and Contact

7.1 Operational-Space Bias Is a Force Vector

For \(D\ddot q+h=\tau\), task \(y=h_y(q)\), and full-row-rank task Jacobian \(J\), suppose the input is implemented as \(\tau=J^\top F+\tau_0\) with \(JD^{-1}\tau_0=0\). Then \[ \Lambda=(JD^{-1}J^\top)^{-1},\qquad \Lambda\ddot y+\eta=F,\qquad \eta=\Lambda(JD^{-1}h-\dot Jv). \] The bias \(\eta\) is a force vector, not a matrix that should be multiplied by task velocity again. Its state dependence generally cannot be expressed using only \(y,\dot y\) in a redundant arm. Rank loss, constraint contacts, torque limits, and imperfect cancellation require separate treatment.

With exact model cancellation, choose \[ F=\eta+\Lambda\left(\ddot y_d-2a\dot e-a^2e\right),\qquad e=y-y_d, \tag{19}\] to obtain \(\ddot e+2a\dot e+a^2e=0\). This is critically damped. Its solutions contain \((c_1+c_2t)e^{-at}\), so a uniform exponential bound with exactly rate \(a\) and fixed prefactor does not hold for arbitrary initial error. Any strictly smaller rate can be bounded with a suitable prefactor. The full error-state metric has dimension twice the task dimension; \(\Lambda\) alone cannot serve that role.

For \(a=1\) and error state \(z=(e,\dot e)\), \[ A_e=\begin{bmatrix}0&1\\-1&-2\end{bmatrix},\qquad P=\begin{bmatrix}3/2&1/2\\1/2&1/2\end{bmatrix} \] satisfies \(A_e^\top P+PA_e=-I\). It gives the valid rate \(1/(2\lambda_{\max}(P))\approx0.292893\ \mathrm{s}^{-1}\) in that metric. In contrast, the simple energy \((e^2+\dot e^2)/2\) has derivative \(-2\dot e^2\), which is zero at \(e\ne0,\dot e=0\) and cannot satisfy a strict instantaneous exponential energy inequality. A conservative certificate can be correct even when it is below the asymptotic spectral rate.

7.2 Saltation and Mode Changes

A contact model needs a guard, pre- and post-contact vector fields, and a reset map. For a smooth non-grazing event \(g(x,t)=0\) with reset \(x^+=R(x^-,t)\), the first-order sensitivity across the event is \[ \Xi=D_xR+ \frac{F^+-D_xR\,F^- -\partial_tR} {D_xg\,F^-+\partial_tg}\,D_xg. \] This is the saltation matrix. It includes event-time shifts as well as reset derivatives; see the saltation review. The denominator must be nonzero, and a consistent coordinate representation and event sequence are required. Grazing, simultaneous contacts, and changes in effective state dimension need additional analysis.

Across an event, a nonexpansion condition is \[ \Xi^\top M^+\Xi\preceq M^-. \] For constant linear mode dynamics, a common flow metric would separately require \[ A_i^\top M+MA_i\preceq-2\lambda M\quad\text{for each mode }i. \tag{20}\] Contracting flows within each mode do not imply this jump condition. The scalar flow \(\dot x=-x\) with a scheduled reset \(x^+=2x^-\) every half second has net multiplier \(2e^{-1/2}>1\). Flow contraction is overwhelmed by repeated expansion at resets. If jumps have bounded expansion rather than nonexpansion, a dwell-time or cumulative flow/jump bound may establish an overall result. A common flow metric alone is insufficient for an impact model.

8 An Executable Certificate Example

The following example verifies the double-integrator identity and metric rate. It is a complete algebraic certificate for the stated linear model, not a DDP solver or a benchmark of golf performance. The generalized symmetric eigenvalue calculation avoids explicitly constructing \(S^{-1/2}\).

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.diag([16.0, 1.0])
R = np.array([[1.0]])
S = solve_continuous_are(A, B, Q, R)
K = np.linalg.solve(R, B.T @ S)
A_closed = A - B @ K
W = Q + K.T @ R @ K
residual = A_closed.T @ S + S @ A_closed + W
rate = float(eigvalsh(W, S)[0] / 2)
assert np.linalg.eigvalsh(S)[0] > 0
assert np.linalg.norm(residual, ord=2) < 1e-10
assert np.isclose(rate, 1.27639320225)
print("Gain:", np.round(K, 6))
print("Metric rate (1/s):", round(rate, 6))

The gain is \([4,3]\) and the metric rate is \(1.276393\) inverse seconds. For a nonlinear model, evaluate the actual policy derivative and the actual metric derivative. Assigning a Riccati derivative to an arbitrary guessed matrix at one state does not establish that a positive metric solves the required differential equation over a trajectory or region.

A numerical study comparing optimizers must publish model parameters, coordinate conventions, integrator and step size, input bounds, costs, initialization, seeds, convergence criteria, metric normalization, feasibility tolerances, and code revisions. Report failed runs and distinguish sampled error fits from worst-direction contraction bounds. Timing requires hardware and solver details. This article supplies no measured success-rate, basin-volume, robot-tracking, or runtime comparison; such results require a reproducible experiment package. An impulse is measured in newton-seconds, whereas a force in newtons requires its duration and profile to define an impulse.

9 A Golf Investigation That Connects the Mechanisms

Start with the delivered outcome: clubhead velocity, face orientation, contact point, and relevant ball-launch variables at the model’s own impact event. A state fixed at a common clock time is a different comparison from a state evaluated when the club reaches the ball. Event-time sensitivity becomes important if perturbed swings arrive earlier or later.

Specify the mechanical state and its history variables. Posture and grip geometry change load transmission. Inertia couples accelerations among links. Earlier muscle work and elastic storage affect the available kinetic energy. Ground reactions can redirect momentum while doing zero work at an ideal stationary contact. Activation and tissue dynamics set feasible corrections and impedance. A same-state affine decomposition describes an acceleration budget; a forward intervention follows a different state and recomputes those mechanisms.

For a candidate feedback policy, propagate the full variational system and contact sensitivities. With state-transition matrix \(\Phi(T,t)\) and a smooth fixed-time outcome \(y_T=H(x(T))\), the first-order output perturbation is \[ \delta y_T=H_x\Phi(T,t)\delta x(t) +\int_t^T H_x\Phi(T,s)B(s)\delta u(s)\,ds. \] For an impact-defined outcome, replace this fixed-time expression by the appropriate event map and include its timing derivatives. The output map, remaining time, disturbance direction, feasible control set, and state uncertainty together determine which errors matter. A metric can express a useful bound, but a large full-state contraction rate is not the only route to small impact variability.

Compare candidate mechanisms using matched outcomes and declared alternatives. A configuration may transmit torque effectively while amplifying a face-orientation error. Co-contraction may reduce one perturbation while increasing energetic cost or changing another response. A low-torque continuation may retain substantial speed without preserving the required contact state. A golfer may allow harmless joint variation while correcting a smaller task-relevant component. These possibilities motivate measurements; none follows solely from an algebraic decomposition or a fitted exponential.

To compare golfers, use the same outcome definition, model class, units, uncertainty treatment, feasible input assumptions, and metric normalization. A computed Riccati matrix depends on chosen cost weights and policy; it is not a directly measured neural strategy. Do not infer that skilled golfers maximize contraction, or that a particular local rate predicts skill, without a preregistered comparison and evidence against alternative explanations. Useful tests include controlled perturbations, repeated trials, parameter sensitivity, grip/contact alternatives, and out-of-sample prediction of delivered impact changes.

10 Exercises and Checks

1. Derive the metric derivative and reproduce the matrix-order counterexample.
2. Compute the alternating discrete product and compare it with its individual spectra.
3. Solve the finite-horizon scalar example and explain why decreasing remaining cost does not imply decreasing state error.
4. Verify the double-integrator ARE, spectral rate, and metric rate separately.
5. Compare the endpoint pullback with the physical two-link mass matrix at straight and bent configurations.
6. Verify the critically damped error-state certificate and identify why its simple energy proof is not strict.
7. Compute the scalar noise variance for common and independent noise realizations.
8. Combine flow and reset multipliers for several event intervals, then include event-time sensitivity in a non-grazing guard example.
9. Define a golf outcome map and identify one state direction that is weakly observed by that outcome. Explain what additional measurements would distinguish low output sensitivity from active error correction.

References and Scope

The contraction definitions follow Lohmiller and Slotine, On Contraction Analysis for Non-linear Systems, Automatica 34(6), 683–696 (1998). The relevant Manchester and Slotine paper is Control Contraction Metrics: Convex and Intrinsic Criteria for Nonlinear Feedback Design, with the amendment linked above. Riccati and DDP identities can be checked against the linked MIT chapters; the numerical counterexamples and certificate code make the distinctions independently reproducible. The saltation review supplies the event-sensitivity framework.

These are mathematical and computational sources. They do not validate a particular golf intervention, neural objective, or physiological efficiency claim. The contribution here is a consistent chain from model assumptions through differential equations and feasible controls to testable outcome predictions.

License: This article and accompanying example code are released under CC BY 4.0.

Cite this page