Full Reference Manuscript: Exact Instantaneous Superposition in Nonlinear Dynamics

Dieter Olson

Full Reference Manuscript: Exact Instantaneous Superposition in Nonlinear Dynamics

Reference Derivations for Tangent Space Linearization and Integral Accumulation

tangent-space

Reference manuscript for the tangent-space series. Use the compact seven-part series as the canonical reading path.

Author

Dieter Olson

Published

January 1, 2026

Reference Status

This is the single-file full reference manuscript for the tangent-space material. It is retained for derivations, appendices, and historical traceability. The canonical rendered reading path is the compact seven-part Tangent-Space Series, which now carries the contraction, hybrid, and residual-aware topics in shorter source-backed form.

Thesis Abstract

A golf swing couples changing posture, segment momentum, shaft deformation, applied torque and contact. Several useful linear relations survive inside this nonlinear motion, but they answer different questions. A differentiable model has a linear derivative at a specified state and input. A control-affine model additionally has an exact finite input-to-rate increment relation at a fixed state. The first variation of a trajectory satisfies a linear differential equation along a common reference motion. None of these facts makes independently evolving swings additive.

Part I develops these distinctions, their geometric meaning, and the regularity needed for error bounds. Part II derives how an entire input perturbation is transported into an endpoint change and how an adjoint weights it by a chosen outcome. Part III develops optimal-control methods subject to their own assumptions. The purpose is a technically usable account of local modeling, with explicit limits on finite predictions and physical interpretation.

A Note on “Exactness”

The derivative and first-variation equation are exact mathematical objects. Replacing a finite nonlinear response by its first-order expansion is an approximation. Continuous differentiability gives a little-o remainder; a quadratic bound needs additional regularity. Exact fixed-state force addition is a separate mechanical statement about a specified input map. Keep these meanings distinct throughout the canonical rendered reading path.

Part I: Geometric Foundations and Local Dynamics

1 Introduction and Motivation

1.1 The Central Thesis

Suppose a model predicts the body’s and club’s future state from initial conditions and an applied-torque history. Three experiments ask three different questions:

At the same posture, velocity and shaft state, how does an extra torque change acceleration now?
Along a specified reference swing, how does a small perturbation change the eventual clubhead state to first order?
If that torque history is actually changed by a finite amount, what new motion and impact follow?

The first can have an exact affine answer. The second is a derivative of the model’s solution map. The third generally requires a new nonlinear integration, including any changed event time. Their relationship explains the usefulness of linear methods without extending their conclusions beyond their hypotheses.

1.2 Central Claim

For a differentiable coordinate model, the map from state/input perturbation directions to the first variation of state rate is linear. If the model is also affine in the declared input, finite input increments produce exactly linear rate increments while state, time and mode remain fixed. Along a common reference, first variations propagate linearly. These are properties of the model and its derivatives; they do not establish biological intent, individual muscle recruitment or optimal swing technique.

1.3 Scope: No Trajectory Superposition

For a homogeneous linear system, linear combinations of solutions are solutions. For a linear system with a fixed nonzero forcing, differences about a common particular solution obey the homogeneous equation. Even in the linear case, adding two complete forced solutions would double the common forcing.

For a nonlinear system, the finite input-to-trajectory map generally fails to be additive. This can happen on ordinary Euclidean space with its usual flat geometry. Curvature of the state manifold is neither required nor a universal measure of the error.

1.4 Terminology: “Tangent Hyperplane” as Affine Subspace

At fixed state $x$, a control-affine model has rates $f_0(x)+G(x)u$. For unrestricted $u\in\mathbb R^m$, this image is the affine subspace $f_0(x)+\operatorname{im}G(x)$ of $T_x\mathcal M$. It is a hyperplane only when its codimension is one. With bounded admissible inputs $\mathcal U(x)$, the actual instantaneous rate set is $f_0(x)+G(x)\mathcal U(x)$, usually only part of that subspace. This set is not the finite-time reachable set.

1.5 Scope and Prerequisites

We use a finite-dimensional state chart and one smooth contact/constraint mode. Assume a unique reference solution exists on the stated finite interval, and all compared solutions remain in the model domain. For the initial derivations, take $f$ continuous in time and continuously differentiable in state and input on a neighborhood of those solutions. Piecewise continuous prescribed inputs can be handled on each smooth time segment. More general measurable inputs require corresponding uniform regularity and integrability hypotheses.

State and input components should have declared units or nondimensional scales before norms are compared. A configuration $q$ is not the complete dynamic state: velocity, shaft modes and retained actuator states matter. At impacts, saturation switches, friction transitions or changing constraint rank, ordinary smooth derivatives may cease to describe the experiment. A smooth equality-constrained system can be reduced to a regular chart; an active inequality or contact transition needs additional analysis.

2 Tangent Spaces and Exact Differentiation

2.1 Nonlinear Dynamics and Smoothness

In local coordinates, let

\[ \dot x=f(x,u,t),\qquad x\in\mathbb R^n,\quad u\in\mathbb R^m. \tag{1}\]

Fix $(\bar x,\bar u,t)$ and combine the perturbations into $h=(h_x,h_u)$. Fréchet differentiability means there is a linear map $D f$ such that

\[ f(\bar x+h_x,\bar u+h_u,t) =f(\bar x,\bar u,t)+D f(\bar x,\bar u,t)[h]+r(h), \tag{2}\]

\[ \lim_{\|h\|\to0}\frac{\|r(h)\|}{\|h\|}=0. \tag{3}\]

The limit controls arbitrary directions approaching zero in the chosen norm. It does not say $r(h)=0$ at a nonzero perturbation. Differentiability at a point does not require a continuous derivative in a whole neighborhood; $C^1$ is a useful stronger hypothesis that gives uniform control on compact neighborhoods and supports the flow results used here.

2.2 Continuous Differentiability Does Not Give a Quadratic Bound

The scalar function $F(z)=|z|^{3/2}$ is $C^1$ near zero and $F'(0)=0$. Its remainder is $|h|^{3/2}$. Dividing by $|h|$ gives a quantity tending to zero, as required for differentiability. Dividing by $h^2$ gives $|h|^{-1/2}$, which is unbounded. Thus a blanket $O(\|h\|^2)$ claim under only $C^1$ assumptions is false.

A locally Lipschitz derivative, often denoted $C^{1,1}$, is sufficient for a quadratic bound. If $\|Df(z+s h)-Df(z)\|\leq L_D s\|h\|$ along the segment, then

\[ r(h)=\int_0^1\bigl(Df(z+s h)-Df(z)\bigr)h\,ds, \qquad \|r(h)\|\leq\tfrac12 L_D\|h\|^2. \]

A continuous bounded second derivative on the relevant neighborhood also supplies this bound. The neighborhood, norm, constant and permissible amplitudes are part of the statement, not optional details.

2.3 Directional Derivatives

For a Fréchet-differentiable map, each directional derivative agrees with the same linear map:

\[ \left.\frac{d}{d\epsilon}f(z+\epsilon h,t)\right|_{\epsilon=0} =Df(z,t)[h]. \tag{4}\]

Existence of derivatives along every line does not by itself establish this property. Define $F(a,b)=a^3/(a^2+b^2)$ away from the origin and $F(0,0)=0$. It is continuous at the origin, and its derivative along direction $(a,b)$ is $a^3/(a^2+b^2)$. The directions $(1,0)$ and $(0,1)$ give 1 and 0, while their sum gives $1/2$. These directional derivatives do not constitute a linear derivative.

2.4 The Jacobian Matrices

Under differentiability, the coordinate derivative has blocks

\[ A(t)=f_x(\bar x(t),\bar u(t),t),\qquad A_{ij}=\frac{\partial f_i}{\partial x_j}, \tag{5}\]

\[ B(t)=f_u(\bar x(t),\bar u(t),t),\qquad B_{ij}=\frac{\partial f_i}{\partial u_j}. \tag{6}\]

The exact first variation of the rate is $A h_x+B h_u$. Both matrices use the same reference state, input and time. A derivative evaluated on one counterfactual branch cannot silently be substituted for that of another.

2.5 A Tangent Space Is Not a Derivative

The tangent space $T_x\mathcal M$ is the vector space of velocities of curves through $x$ on a smooth manifold. It exists independently of any particular dynamic law. A smooth map $F:\mathcal M\to\mathcal N$ has a derivative $T_xF:T_x\mathcal M\to T_{F(x)}\mathcal N$. The space and the map are different objects.

A vector field assigns $f(x)\in T_x\mathcal M$. Geometrically it is a section $\mathcal M\to T\mathcal M$; its derivative therefore maps into the tangent space of the tangent bundle, not automatically back into $T_x\mathcal M$. Local coordinates or an explicitly chosen connection give the familiar component formulas. A state-rate Jacobian along a moving reference is a valid coordinate representation of the variational dynamics; it is not an intrinsic endomorphism independent of that choice.

For a family of solutions $x_\epsilon(t)$ with prescribed $u_\epsilon=\bar u+\epsilon v$ and $x_\epsilon(t_0)=\bar x_0+\epsilon\eta_0$, define

\[ \eta(t)=\left.\frac{\partial x_\epsilon(t)}{\partial\epsilon}\right|_0, \qquad \dot\eta=A(t)\eta+B(t)v(t),\quad\eta(t_0)=\eta_0. \tag{7}\]

This equation follows by differentiating the dynamic equation with respect to the family parameter. $\eta$ is a derivative with respect to amplitude, not the full finite displacement $x_\epsilon-\bar x$. Under the stated regularity on a fixed common interval, the latter is $\epsilon\eta+o(\epsilon)$.

2.6 Coordinate Changes Along a Moving Motion

For a smooth invertible coordinate change $y=\psi(x)$, write $T(t)=D\psi(\bar x(t))$. Variations satisfy $\eta_y=T\eta_x$. Differentiating this identity gives

\[ A_y=(\dot T+T A_x)T^{-1},\qquad B_y=T B_x. \]

The additional $\dot T$ term is essential away from an equilibrium. For $\dot x=1$ on $x>0$ and $y=x^2$, $A_x=0$, while $\dot y=2\sqrt y$ has $A_y=1/x$. A constant $\eta_x$ corresponds to $\eta_y=2x\eta_x$, which grows as $x$ grows. This is the same physical separation represented in different units and coordinates. At an equilibrium of an autonomous system, $\dot T=0$ and the Jacobians are similar, so their eigenvalues agree.

This matters when comparing a joint-space perturbation with a clubhead-space perturbation. A growing coordinate norm need not mean that the physically chosen task error grows at the same rate.

2.7 Example: The Simple Pendulum

For a point mass $m$ at length $\ell$, torque $u$, and angle $\theta$ measured from downward vertical,

\[ m\ell^2\ddot\theta+mg\ell\sin\theta=u. \tag{8}\]

With $x=(\theta,\omega)$,

\[ f(x,u)=\begin{pmatrix}\omega\\-(g/\ell)\sin\theta+u/(m\ell^2)\end{pmatrix}. \tag{9}\]

The derivative is

\[ D f[h_x,h_u]=A h_x+B h_u, \tag{10}\]

\[ A=\begin{pmatrix}0&1\\-(g/\ell)\cos\bar\theta&0\end{pmatrix}, \qquad B=\begin{pmatrix}0\\1/(m\ell^2)\end{pmatrix}. \tag{11}\]

Along any reference solution,

\[ \dot\eta_\theta=\eta_\omega,\qquad \dot\eta_\omega=-(g/\ell)\cos\bar\theta(t)\eta_\theta +v/(m\ell^2). \tag{12}\]

At the unforced downward equilibrium, this becomes

\[ \ddot\eta_\theta+(g/\ell)\eta_\theta=0, \tag{13}\]

\[ \lambda_{1,2}=\pm i\sqrt{g/\ell}. \tag{14}\]

Pure imaginary eigenvalues do not prove stability of the nonlinear equilibrium by the indirect method. Here a separate energy argument does: $E=\tfrac12m\ell^2\omega^2+mg\ell(1-\cos\theta)$ is conserved and has a strict local minimum at $(0,0)$. Sufficiently small sublevel sets remain in a small neighborhood. The equilibrium is locally Lyapunov stable, but nearby positive-energy motions cannot converge to it without dissipation.

For comparison, the planar systems $\dot z=Jz\pm\|z\|^2z$, with $J$ a unit rotation generator, have the same imaginary eigenvalues at the origin. Their radii satisfy $\dot r=\pm r^3$: one decays and the other grows. The derivative alone leaves that distinction unresolved.

The pendulum’s finite-angle discrepancy can be seen from

\[ \sin\theta=\theta-\theta^3/6+O(\theta^5). \tag{15}\]

At this particular equilibrium the first omitted state term is cubic by symmetry. Generic nonlinear models need not have that cancellation.

2.8 Exact Superposition of Input Increments

Linearity of $Df$ gives additivity of derivative directions for any differentiable model. Exact finite input additivity needs the stronger structural condition

\[ f(x,u,t)=f_0(x,t)+G(x,t)u. \]

At fixed $(x,t)$, with one nominal input $\bar u$, define $\Delta_i=f(x,\bar u+v_i,t)-f(x,\bar u,t)$. Then the combined increment is exactly $\Delta_{12}=G(v_1+v_2)=\Delta_1+\Delta_2$. The baseline rate appears once. This statement can hold for finite increments without making the state dependence linear.

For regular mechanical coordinates, $M(q)\ddot q+h(q,\dot q,t)=\mathcal B(q)u$ has this property when $M$ is nonsingular and neither $M$ nor $h$ depends on the declared input. Applied forces are covectors; $M^{-1}\mathcal B$ maps them into acceleration increments. Equality-constraint reactions can preserve affinity after a regular linear solve. A fixed contact mode alone does not guarantee it: a nonlinear actuator law, such as torque proportional to the square of a command, must be inspected in the actual command variable.

If the input set has bounds, the algebraic identity applies whenever the compared inputs are defined; all inputs must also be admissible for a feasible experiment. Bounds do not turn an arbitrary nonlinear actuator into a piecewise linear one.

2.9 Counterfactual Decomposition

Where the zero-input branch is defined,

\[ f(x,u,t)=f(x,0,t)+G(x,t)u. \tag{16}\]

This is a decomposition evaluated at one state. Along an actual motion, integrating the two vector-valued terms in a common coordinate chart reconstructs that motion’s coordinate increment. The integral of $f(x(t),0,t)$ generally does not equal a forward zero-input motion, which evaluates $f$ at its own evolving state. A coordinate increment, a scalar speed, work and a changed impact outcome are different quantities.

3 Time Evolution and Moving Tangent Spaces

3.1 Why Global Superposition Fails

Along a reference motion the coefficients are

\[ A(t)=f_x(\bar x(t),\bar u(t),t),\qquad B(t)=f_u(\bar x(t),\bar u(t),t). \tag{17}\]

Changing the finite input generally changes the entire state history and hence these coefficients. Integrating a derivative on the old history predicts the first variation of the new history, not its complete finite change.

3.2 Flow Transport Is Not Parallel Transport

On $\mathbb R^n$, tangent vectors at different points have a canonical identification as ordinary vectors. Adding them is mathematically meaningful. Nevertheless, for $\dot x=x^2$ on the flat real line,

\[ \varphi_T(x_0)=\frac{x_0}{1-Tx_0},\qquad D\varphi_T(x_0)=\frac{1}{(1-Tx_0)^2}, \]

on an interval before the denominator vanishes. The flow is nonlinear although the space has zero intrinsic curvature. Its derivative transports a perturbation according to the dynamic law. Euclidean parallel transport instead leaves vector components constant; it is a different operation.

More generally, $T\varphi_{t,s}$ transports initial tangent perturbations to the current state along a fixed prescribed input history. It does not require introducing a metric connection. If a connection is used to compare vectors or define norms, state the choice and do not equate its parallel transport with the dynamic solution derivative.

3.4 Variational Equation

The precise local statement is therefore

\[ \dot\eta=A(t)\eta+B(t)v(t),\qquad x_\epsilon(t)=\bar x(t)+\epsilon\eta(t)+o(\epsilon), \tag{20}\]

uniformly on a suitable fixed finite interval. A stronger quadratic remainder requires correspondingly stronger regularity and neighborhood bounds. Verify integration convergence and perturbation-amplitude convergence separately. This is the bridge from exact differentiation to a defensible finite prediction.

The standard flow and equilibrium results are discussed in Hunter’s ODE notes, Sections 1.3–1.4. The distinction between directional and norm-controlled functional variation is developed in Liberzon, Section 1.3.2. The counterexamples above make the specific limitations used here explicit.

4 Local Optimal Control: LQR and Lyapunov Stability

4.1 Connection to Linear Quadratic Regulation

A local derivative model supports controller design when paired with a stated objective and a separate stability argument.

4.1.1 The LQR Problem

Linear Quadratic Regulation (LQR) solves the following problem:

Given:

Linear system: $\dot{x} = Ax + Bu$
Quadratic cost: $J = \int_0^\infty (x^T Q x + u^T R u) dt$
Matrices $Q \succeq 0$ (positive semidefinite), $R \succ 0$ (positive definite)

Find: Optimal feedback control law minimizing $J$

Solution: Linear state feedback $u = -Kx$ where gain $K$ is computed by solving the algebraic Riccati equation:

\[ A^T P + P A - P B R^{-1} B^T P + Q = 0 \tag{21}\]

with $K = R^{-1} B^T P$. Standard sufficient assumptions for the stabilizing infinite-horizon solution are stabilizability of $(A,B)$ and detectability of $(Q^{1/2},A)$. Select the stabilizing positive-semidefinite solution; the matrix equation alone does not guarantee that an arbitrary root defines an optimal stabilizer.

4.1.2 What a Nonlinear LQR Design Establishes

At an equilibrium, a differentiable nonlinear system has a linear first variation. If the resulting closed-loop Jacobian is Hurwitz, the smooth nonlinear closed loop is locally asymptotically stable. This is a local result with a model-dependent neighborhood. It does not establish that the linear feedback globally minimizes the nonlinear objective.

A $C^2$ cost has a second-order Taylor model. Its linear term vanishes at an unconstrained stationary minimum, but the Hessian may be singular or zero: $\ell(x)=x^4$ has no nonzero quadratic leading term. Along a general nominal trajectory, retain cost gradients and state-input cross derivatives rather than assuming each stage is at a minimum.

Along a dynamically feasible moving reference, time-varying LQR propagates a value approximation backward from the terminal cost. Solving independent equilibrium Riccati equations at each instant is a different design. The finite-horizon recursion, nonlinear rollout, input constraints, and stability argument must be checked separately. See the reviewed optimal-control chapter.

4.1.3 Gain Scheduling Requires Its Own Check

Interpolating locally stabilizing gains does not automatically stabilize the scheduled nonlinear system. For a scheduled Lyapunov metric $P(\rho)$, its derivative contains $\sum_i P_{\rho_i}\dot\rho_i$ as well as the closed-loop Jacobian terms. A common Lyapunov function, appropriate scheduling-rate bounds, or a direct trajectory analysis can supply a certificate. Frozen-point eigenvalues alone do not.

4.2 Lyapunov Stability Analysis

4.2.1 Lyapunov’s Indirect Method

One of the most powerful tools for analyzing nonlinear stability is Lyapunov’s indirect method:

Theorem: Lyapunov’s Indirect Method

For $f$ continuously differentiable on a neighborhood, consider $\dot{x} = f(x)$ with equilibrium $x_e$ (i.e., $f(x_e) = 0$). Let $A = \frac{\partial f}{\partial x}(x_e)$ be the Jacobian at equilibrium.

If all eigenvalues of $A$ have strictly negative real parts (Re$(\lambda_i) < 0$), then $x_e$ is locally asymptotically stable
If any eigenvalue of $A$ has strictly positive real part (Re$(\lambda_j) > 0$), then $x_e$ is unstable
If all eigenvalues have non-positive real parts with some on imaginary axis, stability is inconclusive (higher-order analysis required)

4.2.2 Why This Works

For the stable case, write the displacement dynamics as $\dot z=Az+r(z)$, with $r(z)=o(\|z\|)$. If $A$ is Hurwitz, choose $Q\succ0$ and solve $A^\top P+PA=-Q$ for $P\succ0$. Then $V=z^\top Pz$ satisfies

\[ \dot V=-z^\top Qz+2z^\top P r(z). \]

For a sufficiently small neighborhood, the remainder obeys $2\|P\|\|r(z)\|\leq\tfrac12\lambda_{\min}(Q)\|z\|$. Consequently $\dot V\leq-\tfrac12\lambda_{\min}(Q)\|z\|^2$ there. A sufficiently small invariant sublevel set supplies the local stability domain. The proof uses a spectral margin and a bound on the nonlinear remainder; exact differentiation alone is not the certificate.

The Pendulum Example Revisited

Recall the pendulum at $x_e = [0,0]^T$ has:

\[ A = \begin{bmatrix} 0 & 1 \\ -g/\ell & 0 \end{bmatrix}, \quad \lambda = \pm i\sqrt{g/\ell} \]

The eigenvalues are pure imaginary (on the borderline), so Lyapunov’s indirect method is inconclusive. Indeed, the pendulum equilibrium is:

Stable (nearby trajectories stay nearby)
But not asymptotically stable (they don’t converge to equilibrium)
Neutrally stable—consistent with conservative mechanical system

For the inverted pendulum at $x_e = [\pi, 0]^T$, one eigenvalue would be positive, confirming instability.

4.2.3 Design and Validation

The pendulum’s conserved energy, rather than its imaginary eigenvalues alone, establishes the downward equilibrium’s local stability. The inverted equilibrium has a positive eigenvalue and is unstable under the indirect method’s hypotheses.

For a controller, first check the feasible equilibrium or reference and the controllable directions. Then design for the stated cost, derive a closed-loop neighborhood bound and compare its predicted outputs against nonlinear rollouts. Gain scheduling and actuator saturation require their own checks. A small region in a chosen norm does not guarantee good performance if an output is highly sensitive or the control margin is small.

Large excursions, contact transitions and changing shaft states can motivate nonlinear control analysis, but their presence alone does not prove that a particular linear controller will fail. The relevant question is whether the proposed controller meets the declared task tolerances over the required domain.

5 Residuals, Reference Motion, and Nonlinearity

5.1 An Interaction Requires Four Comparable Trajectories

Fix an initial state, model, time interval, and nominal input $\bar u$. Let $x_{00}$ use $\bar u$, $x_{10}$ use $\bar u+\epsilon v_1$, $x_{01}$ use $\bar u+\epsilon v_2$, and $x_{11}$ use $\bar u+\epsilon(v_1+v_2)$. In a common state coordinate chart, define

\[ r(t_1)=x_{11}(t_1)-x_{10}(t_1)-x_{01}(t_1)+x_{00}(t_1). \tag{22}\]

The fourth term is the evolved nominal trajectory, not the initial state. With zero nominal input it is the evolved zero-input trajectory. For the scalar linear system $\dot x=ax+u$, constant $u$, and duration $T$, the endpoint is $e^{aT}x_0+u(e^{aT}-1)/a$ for $a\ne0$. The four-trajectory residual is zero. Replacing its last term by $x_0$ gives $x_0(1-e^{aT})$, a false interaction even in this linear system.

For $\dot x=u^2$ on the Euclidean line and constant perturbations, the same residual is $2\epsilon^2v_1v_2T$. The state space is flat; the input-to-flow map is nonlinear. Ordinary Hessians depend on coordinates and are not intrinsic Riemann curvature tensors. On a manifold, state subtraction requires a declared chart or retraction and a consistent comparison point.

5.2 First Variations and the Mixed Second Variation

For a fixed finite horizon, assume $f(x,u,t)$ is $C^2$ in state and input on a neighborhood containing the relevant trajectories, with sufficient time regularity and bounded input directions. Let

\[ \begin{aligned} \dot\xi_i&=A(t)\xi_i+B(t)v_i,\qquad \xi_i(t_0)=0,\\ A&=f_x(\bar x,\bar u,t),\qquad B=f_u(\bar x,\bar u,t). \end{aligned} \]

These are exact first variations along one nominal trajectory. Their sum obeys the same linear equation. Define $\eta_{12}$ as the mixed second variation of the solution with respect to the two input amplitudes. It satisfies

\[ \dot\eta_{12}=A\eta_{12}+S_{12},\qquad \eta_{12}(t_0)=0, \]

\[ S_{12}=f_{xx}[\xi_1,\xi_2]+f_{xu}[\xi_1,v_2] +f_{xu}[\xi_2,v_1]+f_{uu}[v_1,v_2]. \tag{23}\]

Every derivative is evaluated along the same nominal state and input. For a control-affine system, $f_{uu}=0$, but $f_{xu}$ can be nonzero when the input vector fields depend on state. The interaction is

\[ \begin{aligned} r(t_1)&=\epsilon^2\eta_{12}(t_1)+o(\epsilon^2),\\ \eta_{12}(t_1)&=\int_{t_0}^{t_1}\Phi(t_1,t)S_{12}(t)\,dt. \end{aligned} \tag{24}\]

There is no factor of $1/2$ on this mixed interaction: the two symmetric mixed terms in the Taylor expansion cancel that factor. A single trajectory’s second-order Taylor error is a different quantity.

5.3 Transport Is Part of the Bound

For a chosen consistent norm and fixed interval,

\[ \|r(t_1)\|\leq\epsilon^2\int_{t_0}^{t_1} \|\Phi(t_1,t)\|\,\|S_{12}(t)\|\,dt+o(\epsilon^2). \tag{25}\]

The Jacobian matters through $\Phi$ and the first variations inside $S_{12}$. For $\dot x=ax+u^2$, constant directions about $\bar u=0$ give $r=2\epsilon^2v_1v_2(e^{aT}-1)/a$. A constant Hessian source can therefore be amplified exponentially. A bound proportional to $T$ is valid only when its coefficient controls this transport and the other terms over the entire stated interval; it is not a universal linear growth law.

5.4 What a Small Residual Does Not Establish

A small interaction residual does not bound every prediction error. For instance, two input directions may have no mutual interaction while each has its own nonlinear self-effect. A zero denominator also makes $\|r\|/\|\delta x\|$ unusable. Choose physical scales or a metric, a reference magnitude bounded away from zero, and a task tolerance:

\[ \frac{\|r\|_W}{s_{\mathrm{reference}}}\ll1,\qquad s_{\mathrm{reference}}>0. \tag{26}\]

This is an interaction diagnostic, not a stability or optimality certificate. Separately check finite-model prediction errors, disturbances, discretization error, and active contact modes. A computation should demonstrate convergence under tighter integration tolerances before calling its residual physical.

Part II: Integral Superposition and Transport

6 Fundamentals of Integral Accumulation

6.1 Integration Preserves Linearity

If $g_1,g_2$ are integrable vector-valued functions in one coordinate representation and $a,b$ are fixed scalars, then $\int(ag_1+bg_2)=a\int g_1+b\int g_2$. No dynamic linearity assumption is needed for this algebraic property. The substantive question is whether the integrands being added describe the same experiment.

For first variations along one common reference,

\[ \dot\eta=A(t)\eta+B(t)v(t),\qquad \eta(t_0)=\eta_0. \tag{27}\]

The coefficients are fixed by the reference. The derivative direction $v$ may vary over time. Superposition holds between solutions of this linear variational system when initial variations and input directions are combined consistently.

6.2 The Integration Theorem

Direct integration gives the exact identity

\[ \eta(T)-\eta(t_0) =\int_{t_0}^{T}\bigl(A(t)\eta(t)+B(t)v(t)\bigr)\,dt. \tag{28}\]

The first term includes the unknown evolving variation; omitting it is generally incorrect. Introduce the state-transition matrix $\Phi(t,s)$ by $\partial_t\Phi(t,s)=A(t)\Phi(t,s)$ and $\Phi(s,s)=I$. Variation of constants gives

\[ \eta(T)=\Phi(T,t_0)\eta_0+ \int_{t_0}^{T}\Phi(T,t)B(t)v(t)\,dt. \tag{29}\]

To verify the formula, differentiate its right-hand side with respect to $T$. The integral’s moving upper limit supplies $B(T)v(T)$, and the derivative of each transition matrix supplies $A(T)$ times the complete expression. Its initial value is $\eta_0$. Uniqueness of the linear initial-value problem then identifies it with $\eta$.

6.3 What Superposes: Variations About One Reference

Let $x_{\epsilon,i}$ use input $\bar u+\epsilon v_i$ from the common initial state, and let $x_{00}$ use $\bar u$. The first variation is

\[ \eta_i(T)=\left.\partial_\epsilon x_{\epsilon,i}(T)\right|_0, \qquad x_{\epsilon,i}(T)-x_{00}(T)=\epsilon\eta_i(T)+o(\epsilon). \]

It is not $x_{\epsilon,i}(T)-x(t_0)$. Subtracting the initial state mixes the nominal displacement with the perturbation response. The sum $\eta_1+\eta_2$ is exactly the derivative for $v_1+v_2$; the finite branched responses retain their nonlinear residual.

When an input or initial condition is uncertain, propagate its joint uncertainty through the common derivative. Summing separate output variances is justified only when the required cross-covariances vanish. Linearity of the variational equation does not establish independence of measured joint torques.

7 State Transition Operator

7.1 Solution of the Variational Equation

For zero input variation,

\[ \eta(t)=\Phi(t,s)\eta(s). \tag{30}\]

For finite-dimensional ODEs with continuous coefficients, $\Phi$ is invertible and satisfies $\Phi(t,r)\Phi(r,s)=\Phi(t,s)$ and $\Phi(s,t)=\Phi(t,s)^{-1}$. This is a two-time composition, or cocycle, property. A time-varying system generally depends on both endpoint times, not just their difference.

7.2 Interpretation: Transport Then Accumulate

An input direction applied near time $t$ changes the rate through $B(t)$. The resulting state variation then evolves through $\Phi(T,t)$. A small instantaneous gain can therefore be followed by strong amplification; a large acceleration increment can be followed by cancellation in the selected endpoint output.

Transition matrices at different times need not commute. Take nondimensional matrices

\[ A_1=\begin{pmatrix}0&1\\0&0\end{pmatrix}, \qquad A_2=\begin{pmatrix}0&0\\1&0\end{pmatrix}, \]

on two consecutive unit-duration intervals. Their ordered transition is

\[ e^{A_2}e^{A_1}=\begin{pmatrix}1&1\\1&2\end{pmatrix}, \]

whereas reversing their order gives $\begin{pmatrix}2&1\\1&1\end{pmatrix}$. Neither equals $e^{A_1+A_2}$. Replacing $\Phi$ by an exponential of an integrated Jacobian needs additional commutation conditions.

7.3 Time-Invariant Special Case

For constant $A$,

\[ \Phi(t,s)=e^{A(t-s)}. \tag{31}\]

The LTI solution is

\[ x(T)=e^{A(T-t_0)}x_0+\int_{t_0}^{T}e^{A(T-t)}Bu(t)\,dt. \tag{32}\]

This is the complete motion for the linear plant. The same formula with time-varying reference coefficients gives the derivative response for a nonlinear plant; the identical notation does not make the latter a complete nonlinear solution.

8 Variational Dynamics and Sensitivity Analysis

8.1 Sensitivity to Initial Conditions

Holding the prescribed input history fixed, define

\[ S_x(t)=D_{x_0}x(t;x_0,\bar u). \tag{33}\]

Then

\[ S_x(t)=\Phi(t,t_0), \tag{34}\]

\[ \dot S_x=A(t)S_x,\qquad S_x(t_0)=I. \tag{35}\]

For a differentiable output $y=h(x(T))$, its initial-state derivative is $Dh(\bar x(T))\Phi(T,t_0)$. This task map is essential when the state contains angles, velocities and shaft coordinates while the question concerns clubface orientation or impact position.

8.2 An Input Kernel Is a Density

At fixed terminal time and initial state, write

\[ K(T,t)=\Phi(T,t)B(t). \tag{36}\]

The kernel is a sensitivity density, to be integrated against an input perturbation and time. If input has torque units, its entries have corresponding state-per-torque-per-time units. It is not generally the derivative of the endpoint with respect to changing a finite torque value at one isolated instant.

For bounded measurable input histories, changing the input at a single time leaves its equivalence class almost everywhere unchanged and does not change the integral solution. A short rectangular pulse of height $\epsilon a$ on $[\tau,\tau+h]$ instead produces the first-order endpoint change $\epsilon\int_\tau^{\tau+h}K(T,t)a\,dt$. For continuous $K$ this approaches $\epsilon h K(T,\tau)a$ as $h$ tends to zero.

For the nondimensional system $\dot x=-x+u$, $x(0)=0$, the exact pulse response is

\[ \Delta x(T)=\epsilon a\,e^{-(T-\tau)}(e^h-1). \]

At $T=1$, $\tau=0.3$, $\epsilon a=0.2$ and $h=10^{-4}$, it is about $9.9322\times10^{-6}$, not the height-only value $0.2e^{-0.7}$. A Dirac impulse has a specified area and different units; finite impulsive interventions can leave the small-perturbation regime and need their own dynamic interpretation.

8.3 The Endpoint Derivative Is an Operator

Let $\mathcal E(u)=x(T;x_0,u)$ for a fixed interval and initial state. Under uniform smoothness and neighborhood bounds that make this solution map Fréchet differentiable for bounded input histories, its derivative acts on a whole function:

\[ D\mathcal E(\bar u)[v] =\int_{t_0}^{T}\Phi(T,t)B(t)v(t)\,dt. \tag{37}\]

Writing only $\int\Phi B\,dt$ omits the argument. That matrix is the response to a constant input-offset parameter, not the entire functional derivative.

In the scalar example at $T=1$, $v(t)=1$ gives $1-e^{-1}\approx0.6321$, while $v(t)=t$ gives $e^{-1}\approx0.3679$. Timing changes the endpoint effect even though both directions use the same kernel. This is a precise reason that the sequence of torque and posture changes matters in the swing.

9 Adjoint (Costate) Dynamics

9.1 Backward Sensitivity

For fixed $t_0,T,x_0$, consider

\[ J(u)=\phi(x(T))+\int_{t_0}^{T}L(x(t),u(t),t)\,dt. \tag{38}\]

Here $L$ is the chosen running cost, not necessarily a mechanical Lagrangian. Use column gradients. Along the prescribed nominal history define an adjoint $\lambda$ so that

\[ \begin{aligned} C(t,x)&=\phi(X(T))+\int_t^T L(X(s),\bar u(s),s)\,ds,\\ \lambda(t)&=\left.\nabla_x C(t,x)\right|_{x=\bar x(t)}. \end{aligned} \tag{39}\]

Here $X(s)=x(s;t,x,\bar u)$ starts from state $x$ at time $t$. This is the derivative of the future cost when the future input history is held fixed. It is not automatically the gradient of the optimally minimized value function.

9.2 Adjoint Evolution Equation

The adjoint solves

\[ \dot\lambda=-A^\top\lambda-L_x,\qquad \tag{40}\]

\[ \lambda(T)=\phi_x(\bar x(T)). \tag{41}\]

The first cost variation is

\[ \delta J=\phi_x^\top\eta(T)+ \int_{t_0}^T(L_x^\top\eta+L_u^\top v)\,dt. \]

Using the forward and backward equations, $\frac{d}{dt}(\lambda^\top\eta)=-L_x^\top\eta+\lambda^\top Bv$. Integrating this identity cancels the state-variation terms and yields

\[ DJ(\bar u)[v]=\int_{t_0}^{T} \bigl(L_u+B^\top\lambda\bigr)^\top v\,dt. \tag{42}\]

The quantity $L_u+B^\top\lambda$ is a gradient density for the specified pairing. If the initial state varies, add $\lambda(t_0)^\top\eta_0$. Different inner products on the space of input functions represent the same derivative by different gradients.

For an interior unconstrained stationary input, this density vanishes almost everywhere. With bounds or other constraints, admissible directions are restricted and the appropriate variational inequality or optimal-control conditions replace unrestricted stationarity. Endpoint constraints and variable terminal time also alter boundary conditions.

9.3 The Policy Determines the Sensitivity

If a controller recomputes $u=\kappa(x,t)$ after a state perturbation, the closed-loop state derivative is $A+B\kappa_x$. The running cost becomes $L(x,\kappa(x,t),t)$ and its state gradient includes $\kappa_x^\top L_u$. For $\dot x=x+u$ with $u=-2x$, the closed-loop perturbation decays as $e^{-t}$; freezing the nominal input history instead gives the open-loop sensitivity $e^t$. Both are valid answers to different experiments.

For a golfer, a prescribed torque replay, a model of feedback and a reoptimized motion likewise define different comparisons. An adjoint cannot decide which one represents the player’s response.

9.4 A Reproducible Gradient Check

Use the nondimensional system $\dot x=x^2+u$, $x(0)=0.2$, $T=0.3$, nominal $u=0.4$, $\phi=x^2/2$ and $L=(x^2+0.1u^2)/2$. Then $\dot\lambda=-2x\lambda-x$, $\lambda(T)=x(T)$, and the gradient density is $\lambda+0.1u$. For direction $v(t)=1+t$, compare its integral against $[J(\bar u+\epsilon v)-J(\bar u-\epsilon v)]/(2\epsilon)$ from two independent nonlinear integrations. The supplied regression checks agreement as well as solver accuracy. This verifies the derivative computation for this example; it does not validate a golf model.

9.5 Task and Event Interpretation

An endpoint state sensitivity becomes a clubhead sensitivity only after the appropriate kinematic output derivative is applied. At first contact, the event time usually varies too. For a transversal event $c(x,t)=0$, the time correction has denominator $c_x f+c_t$; a small or zero denominator signals grazing and possible loss of differentiability. Include the event-time derivative, and any impact/reset derivative, before interpreting a pre-impact fixed-time kernel as a ball-strike prediction.

The adjoint weights mechanical effects by a declared objective, while the forward model transports them through changing posture and shaft state. Neither produces a unique allocation of observed performance to muscles, segments or coaching intentions. That inference requires a separately defined intervention and evidence.

10 Physical Interpretations: Impulse and Work

11 Force-Momentum: Impulse Integral

11.1 Newton’s Second Law

For a mechanical Lagrangian $L_m=T-V$ with $p=\partial L_m/\partial\dot q=M(q)\dot q$, the Euler–Lagrange momentum equation can be written:

\[ \frac{dp}{dt} = f(q, \dot{q}) + u \tag{43}\]

Here $u$ is the chosen generalized applied load and $f=\partial L_m/\partial q+Q_{\mathrm{other}}$. The derivative $\partial L_m/\partial q$ includes configuration-dependent kinetic terms as well as potential forces. Canonical generalized momentum need not equal a Cartesian momentum component. For Cartesian translation with constant mass, the equation reduces to Newton’s law with all other external forces included in $f$.

11.2 Impulse Definition

The impulse delivered by a force $u$ over $[t_0, t_1]$ is:

\[ J = \int_{t_0}^{t_1} u(\tau) \, d\tau \tag{44}\]

From Newton’s second law:

\[ \int_{t_0}^{t_1} \frac{dp}{dt} dt = p(t_1) - p(t_0) = \int_{t_0}^{t_1} [f + u] dt \tag{45}\]

If the system starts at rest ($\dot{q}(t_0) = 0$) and the integral of every retained term in $f$ is demonstrably negligible on the stated interval, then approximately:

\[ p(t_1) \approx \int_{t_0}^{t_1} u(\tau) d\tau = J \tag{46}\]

The integral of the complete right-hand side equals the canonical momentum change. The impulse of the selected input alone equals that change only when all remaining integrated terms vanish; in Cartesian coordinates the total external impulse equals the change in linear momentum.

11.3 Superposition of Impulses

If two forces $u_1$ and $u_2$ are applied:

\[ \begin{aligned} J_{12}&=\int_{t_0}^{t_1}[u_1(\tau)+u_2(\tau)]\,d\tau\\ &=\int_{t_0}^{t_1}u_1(\tau)\,d\tau+ \int_{t_0}^{t_1}u_2(\tau)\,d\tau=J_1+J_2. \end{aligned} \tag{47}\]

This is exact—impulses superpose because integration is linear. The resulting momentum change is:

\[ p(t_1) = p(t_0) + J_1 + J_2 + \int_{t_0}^{t_1}f(q(t),\dot q(t))\,dt \tag{48}\]

This is exact bookkeeping on the combined trajectory. The remaining integral is not generally higher order or small; gravity can contribute at first order in elapsed time. If force laws or feedback depend on the state, their actual histories also change across separate counterfactual runs.

Interpretation

The complete momentum balance adds retained load terms in one coordinate representation along the actual motion. Its individual input integrals are signed accounting terms, not independent forward-trajectory predictions.

12 Power-Energy: Work Integral

12.1 Power and Work

For a mechanical system, power delivered by a force is:

\[ P(t) = u(t)^T \dot{q}(t) \tag{49}\]

The work done over $[t_0, t_1]$ is:

\[ W = \int_{t_0}^{t_1} P(\tau) d\tau = \int_{t_0}^{t_1} u(\tau)^T \dot{q}(\tau) d\tau \tag{50}\]

For generalized load $u$ paired with $\dot q$, this is the work of that load. The kinetic-energy balance must include the work $W_{\mathrm{other}}$ of all other forces, including gravity when it is not absorbed into potential energy:

\[ W+W_{\mathrm{other}}=E_k(t_1)-E_k(t_0) \tag{51}\]

Here $E_k=\tfrac12\dot q^TM\dot q$. For a time-independent conservative potential and only dissipative loads $-d$ besides $u$, the mechanical-energy balance is $\Delta(E_k+V)=\int(u-d)^T\dot q\,dt$. Ideal stationary constraint reactions do no work; moving supports can. These assumptions must be stated before equating actuator work with an energy change.

12.2 Superposition of Work

If two forces $u_1$ and $u_2$ are applied:

\[ \begin{aligned} W_{12}&=\int_{t_0}^{t_1}[u_1+u_2]^\top\dot q\,d\tau\\ &=\int_{t_0}^{t_1}u_1^\top\dot q\,d\tau+ \int_{t_0}^{t_1}u_2^\top\dot q\,d\tau=W_1+W_2. \end{aligned} \tag{52}\]

Again, exact additivity due to integral linearity.

The velocity history is the same combined-motion history in both work terms. For a unit mass initially at rest, two constant unit forces over a unit interval produce velocity $2t$, final kinetic energy 2 and work 1 from each force. Running either force alone gives velocity $t$ and work $1/2$. The separate-run works sum to 1, not 2. This interaction occurs even with linear force-to-velocity dynamics because power multiplies force by velocity.

Variation vs. Causation

The force-work terms add along the same combined trajectory, providing a valid signed energy budget when all terms are included. Those terms are not the works obtained by running each input alone, and normalized work fractions do not uniquely identify physiological causation or efficiency. Negative work and an initial energy store can also make naive percentage interpretations misleading.

Work accounting and first-variation transport are distinct linear operations. Neither automatically assigns a unique fraction of an observed outcome to a physical cause.

13 General Principle: Integral Accumulation

Both impulse and work illustrate a general principle:

Integral Accumulation Principle

Along a fixed trajectory $x(t)$, for any quantity defined as an integral of a linear function of inputs:

\[ \Psi = \int_{t_0}^{t_1} g(x(t), u(t), t) \, dt \]

where $g$ is linear in $u$ with $x(t)$ held fixed, the contributions from different inputs superpose exactly:

\[ \Psi(u_1 + u_2) = \Psi(u_1) + \Psi(u_2) \]

Important caveat: This integral identity holds for any such input histories evaluated along the same fixed trajectory $x(t)$; the histories need not be infinitesimal. If $x(t)$ itself depends nonlinearly on $u$ (as it generally does in closed-loop or trajectory optimization settings), then the principle does not hold globally — changing $u$ changes the trajectory, and the two contributions cannot be evaluated independently.

Impulse and work are special cases where $g = u$ and $g = u^T \dot{q}$ respectively. The linearity of integration ensures exact additivity of these quantities, even as the underlying dynamics remain nonlinear.

14 Finite-Horizon Residuals and Error Quantification

14.1 Keep the Baseline and the Perturbation Definition Fixed

The four-trajectory interaction in Equation 22 can also be written using finite deviations $\Delta x_{ij}=x_{ij}-x_{00}$:

\[ r(t_1)=x_{11}(t_1)-x_{10}(t_1)-x_{01}(t_1)+x_{00}(t_1). \tag{53}\]

\[ r=\Delta x_{11}-\Delta x_{10}-\Delta x_{01}. \tag{54}\]

These finite deviations differ from the first variations $\xi_i$, whose equation is linear and whose interaction is exactly zero. The notation must not quietly turn one into the other.

For a smooth solution map on a fixed interval,

\[ \|r\|=O(\epsilon^2). \tag{55}\]

If $\|\Phi(t_1,t)\|\,\|S_{12}(t)\|\leq C_T$ throughout an interval of duration $T$, then

\[ \|r\|\leq\epsilon^2 C_T T+o(\epsilon^2). \tag{56}\]

$C_T$ can grow rapidly with the horizon. Quadratic order does not mean exact quadrupling at arbitrary amplitudes: the mixed coefficient can vanish, higher orders can matter, and contacts can destroy smooth dependence. Convergence in perturbation amplitude should be tested at a fixed horizon.

14.2 Interaction and Single-Trajectory Approximation Error

The relevant Hessians are

\[ H_x=f_{xx},\qquad H_u=f_{uu},\qquad H_{xu}=f_{xu}. \tag{57}\]

The mixed interaction uses $S_{12}$ from Equation 23:

\[ r(t_1)=\epsilon^2\int_{t_0}^{t_1}\Phi(t_1,t)S_{12}(t)\,dt+o(\epsilon^2). \tag{58}\]

By comparison, a single-input expansion $x_\epsilon=\bar x+\epsilon\xi+\tfrac12\epsilon^2\eta+o(\epsilon^2)$ has a self-variation source $f_{xx}[\xi,\xi]+2f_{xu}[\xi,v]+f_{uu}[v,v]$. Both calculations need transport. Neither identifies intrinsic manifold curvature from a measured residual.

14.3 A Task-Based Interpretation for the Swing

Specify whether comparisons use the same elapsed time or an event such as first ball contact. A perturbed trajectory may reach that event at a different time, so a contact-event endpoint map must include event-time sensitivity; smooth-flow expressions alone do not provide it. A changed contact mode requires a hybrid analysis.

In a common coordinate chart, compare residuals using declared scales:

\[ \chi=\|r\|_W/s_{\mathrm{reference}},\qquad s_{\mathrm{reference}}>0. \tag{59}\]

If the chosen perturbation response scale is proportional to $\epsilon$, a bounded mixed response gives a diagnostic of order

\[ \chi=O(\epsilon C_T T). \tag{60}\]

For a prescribed tolerance $\chi_{\max}$, a leading-order estimate is

\[ \epsilon C_T T\lesssim\chi_{\max}, \tag{61}\]

with scales absorbed into $C_T$. It is meaningful only where the derivative bounds and remainder estimates hold. There is no universal threshold such as 0.1 at which a nonlinear controller becomes necessary.

Applied to a swing model, this procedure can quantify how two chosen input interventions interact in club position, face orientation, or speed. It does not recover a unique causal allocation of observed performance to muscles. Geometry, state-dependent leverage, activation, shaft storage, feedback, and impact timing can all alter the endpoint map. Compare those mechanisms through explicit model variants and measured uncertainties.

Part III: Hamiltonian Structure and Optimal Control

15 Variational Principles and Hamiltonian Formulation

15.1 Mechanics and Optimization Ask Different Questions

Hamilton’s stationary-action principle determines conservative mechanical trajectories by varying configurations with fixed endpoints in the action $\int L_m(q,\dot q,t)\,dt$. The stationary action need not be a minimum. Applied nonconservative forces enter through the Lagrange–d’Alembert formulation.

Optimal control instead chooses an admissible input to minimize a specified objective subject to dynamics. Its running cost need not be kinetic minus potential energy, and the mechanical system does not generally minimize it without an appropriately designed input. Here use $\ell$ for that cost and $L_m$ for a mechanical Lagrangian.

For a fixed initial state, the objective is

\[ J(u)=\phi(x(T),T)+\int_{t_0}^T\ell(x,u,t)\,dt. \tag{62}\]

The state equation is $\dot x=f(x,u,t)$, and admissible inputs, endpoint conditions and path constraints must be specified. All gradients below are column vectors.

15.2 A Consistent Minimization Convention

For a normal extremal, define

\[ H(x,u,\lambda,t)=\ell(x,u,t)+\lambda^\top f(x,u,t). \tag{63}\]

The corresponding necessary conditions are

\[ \dot x=H_\lambda=f,\qquad \dot\lambda=-H_x,\qquad u^*(t)\in\underset{u\in\mathcal U}{\operatorname{argmin}}H. \tag{64}\]

The pointwise condition holds almost everywhere under the hypotheses of the relevant optimal-control principle. For an unconstrained interior solution, it implies $H_u=0$. Stationarity is weaker than minimization, and necessary conditions alone do not prove optimality.

Some texts use a Hamiltonian with negative running cost and maximize it. For example, Liberzon, Section 3.4, uses $p^\top f-\ell$. Setting $p=-\lambda$ makes that Hamiltonian the negative of ours. Both conventions are consistent; combining our positive-cost Hamiltonian with that maximization rule reverses the result.

The general multiplier formulation includes a nonnegative cost multiplier $\lambda_0$ in $H=\lambda_0\ell+\lambda^\top f$, with nontrivial multipliers. The normalization $\lambda_0=1$ applies to normal extremals. Abnormal extremals with $\lambda_0=0$ require separate treatment; they are not excluded merely by writing the normal equations.

15.3 Costate as Sensitivity

Along a fixed nominal input history, the adjoint represents the gradient of the future cost defined in Part II. At an optimum, it can also agree with the gradient of a differentiable value function when the associated dynamic-programming conditions hold. Neither interpretation makes a general costate a physical momentum.

For a fixed terminal time and a free terminal state with only terminal cost,

\[ \lambda(T)=\phi_x(x(T),T). \tag{65}\]

Backward integration also needs the nominal state and input histories: a terminal state and an objective alone do not supply all coefficients of the adjoint equation.

15.4 Deriving the Dynamic Residual Correctly

Consider the fixed-time, unconstrained-interior problem

\[ \min_u\ \phi(x(T),T)+\int_{t_0}^T\ell(x,u,t)\,dt \tag{66}\]

subject to

\[ \dot x=f(x,u,t),\qquad x(t_0)=x_0. \tag{67}\]

Introduce a multiplier for the complete dynamic residual:

\[ \widetilde J=\phi(x(T),T)+ \int_{t_0}^T\left[\ell+\lambda^\top(f-\dot x)\right]\,dt. \tag{68}\]

With

\[ H=\ell+\lambda^\top f, \tag{69}\]

this is

\[ \widetilde J=\phi+\int_{t_0}^T(H-\lambda^\top\dot x)\,dt. \tag{70}\]

Vary $x,u,\lambda$ independently, with $\delta x(t_0)=0$. The first variation is

\[ \begin{aligned} \delta\widetilde J={}&\phi_x^\top\delta x(T)\\ &+\int_{t_0}^T\left[ H_x^\top\delta x+H_u^\top\delta u +(f-\dot x)^\top\delta\lambda -\lambda^\top\delta\dot x\right]dt. \end{aligned} \tag{71}\]

The coefficient of $\delta\lambda$ is $f-\dot x$, not just $H_\lambda=f$. Omitting the second term would impose $f=0$ instead of the state equation.

Integration by parts gives

\[ -\int_{t_0}^T\lambda^\top\delta\dot x\,dt =-\left[\lambda^\top\delta x\right]_{t_0}^{T} +\int_{t_0}^T\dot\lambda^\top\delta x\,dt. \tag{72}\]

Consequently,

\[ \begin{aligned} \delta\widetilde J={}&(\phi_x-\lambda(T))^\top\delta x(T)\\ &+\int_{t_0}^T\left[ (H_x+\dot\lambda)^\top\delta x +H_u^\top\delta u +(f-\dot x)^\top\delta\lambda\right]dt. \end{aligned} \tag{73}\]

Vanishing coefficients yield the state equation, costate equation, interior stationarity and free-state terminal condition. A fixed terminal state has $\delta x(T)=0$, so its multiplier is not determined by that free-state condition. Bounded controls restrict the admissible variations; needle variations underpin the stronger pointwise minimum principle.

15.5 Endpoint and Time Conditions

With a regular endpoint equality $\psi(x(T),T)=0$, an endpoint multiplier $\nu$ gives

\[ \lambda(T)=\phi_x+\psi_x^\top\nu. \]

When terminal time is also free and the usual regularity conditions apply,

\[ H(T)+\phi_T+\nu^\top\psi_T=0. \]

These are partial derivatives of the explicit endpoint functions. For an autonomous minimum-time problem with fixed terminal state and no terminal cost, the normal Hamiltonian is zero along an extremal. A fixed terminal time does not supply this zero-Hamiltonian condition.

15.6 Bounded Inputs and Bang-Bang Control

For our minimization convention,

\[ u^*(t)\in\underset{u\in\mathcal U}{\operatorname{argmin}} H(x^*,u,\lambda^*,t). \tag{74}\]

If $H$ is affine in a scalar $u\in[-1,1]$, a nonzero coefficient selects an endpoint. If the coefficient vanishes over an interval, the minimum condition alone does not select the input; a singular arc may require higher-order conditions. Thus bounded inputs do not universally imply a specified bang-bang sequence.

15.7 Minimum Squared Input Is Not Mechanical Energy

For the LTI system $\dot x=Ax+Bu$, fixed endpoints and $R\succ0$, consider

\[ J=\tfrac12\int_{t_0}^T u^\top Ru\,dt. \tag{75}\]

This is a weighted squared-input objective. Calling it energy requires a separate physical model relating the input to power. The Hamiltonian and stationary input are

\[ H=\tfrac12u^\top Ru+\lambda^\top(Ax+Bu), \tag{76}\]

\[ u^*=-R^{-1}B^\top\lambda, \tag{77}\]

\[ \dot\lambda=-A^\top\lambda. \tag{78}\]

The state/costate boundary-value system is

\[ \frac{d}{dt}\begin{pmatrix}x\\\lambda\end{pmatrix} =\begin{pmatrix}A&-BR^{-1}B^\top\\0&-A^\top\end{pmatrix} \begin{pmatrix}x\\\lambda\end{pmatrix}. \tag{79}\]

Let $d=x_T-e^{A(T-t_0)}x_0$ and $W=\int_{t_0}^T e^{A(T-t)}BR^{-1}B^\top e^{A^\top(T-t)}\,dt$. If $W$ is nonsingular, the optimal input is $u^*(t)=R^{-1}B^\top e^{A^\top(T-t)}W^{-1}d$ and $J^*=\tfrac12d^\top W^{-1}d$. For singular $W$, feasibility requires $d$ in its range; the corresponding minimum uses the pseudoinverse on that reachable subspace.

For the nondimensional double integrator from rest to rest over displacement $D$ in fixed time $T$, this gives $u(t)=6D/T^2-12Dt/T^3$ and $J=6D^2/T^3$ when $R=1$. The checked example $D=3,T=2$ has zero final velocity, displacement 3 and cost 6.75. If input bounds exclude that history, it is not the constrained optimum.

15.8 A Complete Minimum-Time Example

For $\dot q=v$, $\dot v=u$, $|u|\leq1$, transfer $(q,v)=(0,0)$ to $(1,0)$ in minimum nondimensional time:

\[ J=T=\int_0^T1\,dt, \tag{80}\]

\[ H=1+\lambda_qv+\lambda_vu. \tag{81}\]

The minimizing rule is

\[ u^*=-\operatorname{sign}(\lambda_v),\qquad\lambda_v\ne0. \tag{82}\]

The solution uses $u=1$ for $0<t<1$ and $u=-1$ for $1<t<2$. Its velocity is $t$ then $2-t$, so the traveled displacement is 1. For any rest-to-rest history of duration $T$, $v(t)\leq\min(t,T-t)$ and hence $q(T)\leq T^2/4$. The achieved duration 2 therefore meets a global lower bound for this specific problem.

Choose $\lambda_q=-1$, $\lambda_v=t-1$. Then

\[ H=1-v+(t-1)u=0 \tag{83}\]

on both arcs, verifying the normal free-time condition. The input value at the isolated switching instant does not change the trajectory. This example’s single switch is not a theorem about arbitrary nonlinear plants.

15.9 LQR as a Special Boundary-Value Problem

For

\[ J=\tfrac12x(T)^\top Q_f x(T) +\tfrac12\int_{t_0}^T(x^\top Qx+u^\top Ru)\,dt, \tag{84}\]

\[ \dot x=Ax+Bu, \tag{85}\]

the Hamiltonian is

\[ H=\tfrac12(x^\top Qx+u^\top Ru)+\lambda^\top(Ax+Bu). \tag{86}\]

With $Q,Q_f\succeq0$ and $R\succ0$,

\[ u^*=-R^{-1}B^\top\lambda, \tag{87}\]

\[ \dot\lambda=-Qx-A^\top\lambda. \tag{88}\]

Substituting $\lambda=P(t)x$ and matching coefficients yields

\[ -\dot P=Q+A^\top P+PA-PBR^{-1}B^\top P, \qquad P(T)=Q_f, \tag{89}\]

\[ u^*=-R^{-1}B^\top P(t)x. \tag{90}\]

This solves the stated finite-horizon linear-quadratic problem. For a nonlinear plant it supplies a local approximation or controller design, whose finite performance and stability require separate verification.

15.10 Box Constraints Require the Correct Metric

For $u_{\min}\leq u\leq u_{\max}$, use inequalities $u_{\min}-u\leq0$ and $u-u_{\max}\leq0$. The stationarity condition is

\[ H_u-\mu_{\min}+\mu_{\max}=0, \tag{91}\]

with nonnegative multipliers and componentwise complementarity,

\[ \mu_{\min}\odot(u-u_{\min})=0,\qquad \mu_{\max}\odot(u_{\max}-u)=0. \tag{92}\]

For a convex quadratic Hamiltonian, minimization is projection of the unconstrained solution in the quadratic weight’s metric. Componentwise clipping solves this projection for a positive diagonal weight and box constraints; it generally fails for a coupled weight.

For example, minimize $\tfrac12(u-a)^\top R(u-a)$ with $a=(2,0)$, $R=\begin{pmatrix}2&1\\1&2\end{pmatrix}$ and $u\in[-1,1]^2$. Clipping gives $(1,0)$, while the actual optimum is $(1,1/2)$. Its gradient is $(-3/2,0)$, balanced by upper-bound multiplier $(3/2,0)$. Changing one saturated component changes the best value of the other through the cross term.

15.11 Path Constraints and Their Limits

For regular mixed state/input inequalities $c(x,u,t)\leq0$, a smooth active-arc multiplier formulation can use

\[ H_{\mathrm{aug}}=H+\mu^\top c, \tag{93}\]

\[ \dot\lambda=-H_x-c_x^\top\mu, \tag{94}\]

together with stationarity, feasibility, nonnegative multipliers and complementarity. Applicable constraint qualifications are required.

Pure state constraints can involve boundary arcs, differentiated tangency conditions and costate jumps or measure multipliers. They are not universally handled by appending an ordinary smooth multiplier to these equations. A numerical transcription must check constraints between nodes as well as at nodes, and a contact reset needs its own model. These distinctions matter whenever a swing model imposes grip closure, ground contact, actuator bounds or a precise impact event.

16 Numerical Implementation and Discrete-Time Systems

16.1 Differentiate the Actual Discrete Transition

Let $F_k(x,u)$ be the transition used by the implementation, with an explicitly specified input hold over $[t_k,t_{k+1}]$. It may be the exact held-input flow or a numerical approximation. A nominal rollout satisfies $\bar x_{k+1}=F_k(\bar x_k,\bar u_k)$. Euler defines

\[ F_k^E(x,u)=x+h_kf(x,u,t_k). \tag{95}\]

The expansion about that nominal transition is

\[ \begin{aligned} x_{k+1}={}&\bar x_{k+1}+A_k^d\delta x_k+B_k^d\delta u_k\\ &+o(\|(\delta x_k,\delta u_k)\|). \end{aligned} \tag{96}\]

Here

\[ A_k^d=D_xF_k,\qquad B_k^d=D_uF_k. \tag{97}\]

For Euler these are $I+h_kf_x$ and $h_kf_u$, respectively. For an RK4 rollout, differentiate its entire transition, including all intermediate evaluations. Using Euler derivatives with an RK4 rollout is an additional approximation, not the exact derivative of the implemented model. An arbitrarily drawn state path has defects $F_k(\bar x_k,\bar u_k)-\bar x_{k+1}$; either include those defects in the optimization or generate the nominal states by rollout.

The first variation of the chosen discrete model obeys exactly

\[ \eta_{k+1}=A_k^d\eta_k+B_k^dv_k. \tag{98}\]

This exact derivative statement does not assert exact finite perturbation prediction. Even a continuous system affine in held input can have a nonlinear finite-step map in that input.

16.2 Discrete Transport and Input Holding

Define

\[ \Phi_{k+1|k}=A_k^d, \tag{99}\]

\[ \Phi_{N|0}=A_{N-1}^d\cdots A_1^dA_0^d. \tag{100}\]

Then $\eta_N=\Phi_{N|0}\eta_0+\sum_{k=0}^{N-1}\Phi_{N|k+1}B_k^dv_k$. The empty product is the identity. The order matters.

For $\dot q=v,\ \dot v=u$ with constant input over duration $h$, the exact matrices are $A^d=\begin{pmatrix}1&h\\0&1\end{pmatrix}$ and $B^d=(h^2/2,h)^\top$. Euler instead gives $(0,h)^\top$ for the input column. The difference is real displacement accumulated while velocity changes; it is not a change in the dimension of a tangent space.

16.3 Integration Order and Regularity

For autonomous dynamics with a held input, classical RK4 uses

\[ \begin{aligned} k_1&=f(x_k,u_k),\\ k_2&=f(x_k+\tfrac h2k_1,u_k),\\ k_3&=f(x_k+\tfrac h2k_2,u_k),\\ k_4&=f(x_k+hk_3,u_k),\\ F_k^{RK4}&=x_k+\tfrac h6(k_1+2k_2+2k_3+k_4). \end{aligned} \tag{101}\]

For explicitly time-dependent dynamics, use the corresponding stage times. Input interpolation must also match the specified hold. On a fixed smooth interval with sufficient bounded derivatives and stable error propagation, Euler has one-step state error $O(h^2)$ and global error $O(h)$; RK4 has one-step error $O(h^5)$ and global error $O(h^4)$. Some references divide the one-step defect by $h$ when defining local truncation error; the convention must be stated.

An input jump or impact can violate the smoothness needed by a high-order formula. Align integration segments with known discontinuities and locate events to a stated tolerance. Increasing formal order does not repair an incorrect contact model.

16.4 Adaptive Integration Is Not a Stiffness Cure

An embedded pair estimates local error from two formulas. A simplified controller for an estimate proportional to $h^5$ is

\[ h_{\rm new}=s\,h_{\rm old} \left(\frac{\mathrm{tol}}{e_k}\right)^{1/5}, \tag{102}\]

with safety factor $s<1$, bounded step changes, component scaling and safeguards for zero error. Production solvers combine relative and absolute tolerances; position and angular velocity need not share one unscaled tolerance.

Dormand–Prince RK45 is an explicit nonstiff method. SciPy’s solver documentation distinguishes explicit RK methods from Radau/BDF options for stiff problems. The requested output times are not necessarily the internal step sequence, and event detection can miss multiple crossings within one step.

For held inputs, backward Euler and implicit trapezoidal transitions satisfy

\[ x_{k+1}=x_k+h f(x_{k+1},u_k,t_{k+1}), \tag{103}\]

\[ \begin{aligned} x_{k+1}=x_k+\tfrac h2[& f(x_k,u_k,t_k)\\ &+f(x_{k+1},u_k,t_{k+1})]. \end{aligned} \tag{104}\]

Their nonlinear equations must be solved sufficiently accurately. Differentiate the converged implicit relation, including its inverse Jacobian, when forming transition sensitivities.

16.5 A Convergence Bound With Explicit Assumptions

Suppose the trajectory remains in a domain where

\[ \|f(x_1,u,t)-f(x_2,u,t)\|\leq L\|x_1-x_2\|, \tag{105}\]

and its second time derivative is bounded by $M$. With exact initial state and uniform Euler step $h$, the recursion $e_{k+1}\leq(1+Lh)e_k+Mh^2/2$ yields

\[ e_N\leq\frac{M}{2L}(e^{L(T-t_0)}-1)h. \tag{106}\]

For $L=0$, use the limit $M(T-t_0)h/2$. This is a finite-horizon propagation argument under the stated domain and smoothness assumptions. It does not need a claim that physical perturbations never grow exponentially, or an appeal to a PDE equivalence theorem.

For example, Euler converges to the physically unstable solution of $\dot x=x$. Halving $h$ asymptotically halves its fixed-time error. For RK4 the global error ratio approaches 16, whereas the one-step error ratio approaches 32. These are different experiments.

16.6 Absolute Stability and Fast Modes

For the scalar test equation $\dot x=\lambda x$, a method has amplification $R(h\lambda)$. Its absolute-stability region satisfies $|R(z)|\leq1$. Euler has $R(z)=1+z$; backward Euler has $R(z)=1/(1-z)$; trapezoidal integration has $R(z)=(1+z/2)/(1-z/2)$.

Consider

\[ \dot x_1=-0.01x_1,\qquad \dot x_2=-1000x_2. \tag{107}\]

Euler requires $h<0.002$ for strict decay of the fast mode. The value $0.002$ is the neutral boundary, not a desirable accuracy choice. Backward Euler and trapezoidal integration are A-stable, but only backward Euler among these two is L-stable. At $z=-100$, their amplifications are $1/101$ and $-49/51$: the latter barely damps an extremely fast physical decay and alternates sign. Implicit does not mean universally stable, and A-stability does not imply accurate resolution of every mode.

16.7 Keep the Error Sources Separate

Under suitable bounds, accumulated integration error can be estimated by

\[ e_N\leq C N h^{p+1}e^{LNh}, \tag{108}\]

so, for fixed $Nh=T$,

\[ e_N=O(h^p). \tag{109}\]

In contrast, for a smooth endpoint map and a finite input perturbation of size $\epsilon$,

\[ \mathcal E_{\rm variation}(T,\epsilon)=O(\epsilon^2) \tag{110}\]

when its second derivative is bounded. This error generally approaches a nonzero value as $h\to0$ at fixed $\epsilon$. Part I’s $\dot x=x^2$ example demonstrates it exactly. Repeated optimization updates can improve a nominal trajectory, but no universal per-iteration $O(h)$ prediction-error theorem follows.

For a swing calculation, report integration convergence, perturbation-amplitude convergence, optimization residuals, constraint/event errors and uncertainty in the mechanical model separately. A tiny integrator tolerance cannot establish that the inferred joint forces or shaft parameters are correct.

16.8 Verify Transition Derivatives

The quantities needed by the discrete optimizer are

\[ A_k^d=D_xF_k,\qquad B_k^d=D_uF_k. \tag{111}\]

A central directional check is

\[ D_xF_k\,v\approx \frac{F_k(x+\epsilon v,u)-F_k(x-\epsilon v,u)}{2\epsilon}. \tag{112}\]

Sweep several scales of $\epsilon$ to reveal truncation and roundoff effects, using dimensionally scaled variables. Automatic differentiation differentiates the executed computational graph to floating-point accuracy; it does not remove model error, nondifferentiable branches, solver tolerance effects or a missing event-time derivative. Compare against an independent calculation before interpreting a gradient as physical evidence.

17 Control Synthesis via DDP, iLQR, and MPC

17.1 One Quadratic Model for DDP and iLQR

Both differential dynamic programming (DDP) and iterative LQR (iLQR) use a quadratic local cost model along a nominal nonlinear rollout. Full DDP also retains second derivatives of the discrete dynamics in the local action-value Hessian; iLQR omits those particular terms. A first-order state-cost expansion is not the defining iLQR approximation.

Use $F_k$ for the actual discrete transition and $L_k$ for the discrete running cost, including the chosen quadrature and timestep. For fixed initial state,

\[ V_k(x)=\min_{u_k,\ldots,u_{N-1}} \left[\phi(x_N)+\sum_{j=k}^{N-1}L_j(x_j,u_j)\right], \tag{113}\]

subject to $x_{j+1}=F_j(x_j,u_j)$. Bellman’s equation is

\[ \begin{aligned} V_k(x)&=\min_u[L_k(x,u)+V_{k+1}(F_k(x,u))],\\ V_N&=\phi. \end{aligned} \tag{114}\]

General nonlinear value functions need not be smooth everywhere. The expansions below describe a smooth local model near a chosen nominal trajectory, not a globally quadratic value function.

17.2 Derivatives of the Local Action Value

At $\bar x_{k+1}=F_k(\bar x_k,\bar u_k)$, write

\[ V_{k+1}(\bar x_{k+1}+\eta) \approx V^+ + (v^+)^\top\eta+\tfrac12\eta^\top P^+\eta. \tag{115}\]

With $A=D_xF_k$, $B=D_uF_k$,

\[ \eta=A\delta x+B\delta u+o(\|(\delta x,\delta u)\|). \tag{116}\]

The running-cost expansion retains all quadratic blocks:

\[ \begin{aligned} L_k(\bar x+\delta x,\bar u+\delta u)\approx{}& L^0+L_x^\top\delta x+L_u^\top\delta u\\ &+\tfrac12\delta x^\top L_{xx}\delta x +\delta x^\top L_{xu}\delta u\\ &+\tfrac12\delta u^\top L_{uu}\delta u. \end{aligned} \tag{117}\]

Define the action value using the same absolute-state transition,

\[ \mathcal Q_k(x,u)=L_k(x,u)+V_{k+1}(F_k(x,u)). \tag{118}\]

Its local approximation is

\[ \begin{aligned} \mathcal Q_k\approx{}&q^0+q_x^\top\delta x+q_u^\top\delta u\\ &+\tfrac12\delta x^\top Q_{xx}\delta x +\delta x^\top Q_{xu}\delta u\\ &+\tfrac12\delta u^\top Q_{uu}\delta u. \end{aligned} \tag{119}\]

The chain rule gives full-DDP derivatives:

\[ \begin{aligned} q_x&=L_x+A^\top v^+,\\ q_u&=L_u+B^\top v^+,\\ Q_{xx}&=L_{xx}+A^\top P^+A+\sum_i v_i^+(F_{k,i})_{xx},\\ Q_{xu}&=L_{xu}+A^\top P^+B+\sum_i v_i^+(F_{k,i})_{xu},\\ Q_{uu}&=L_{uu}+B^\top P^+B+\sum_i v_i^+(F_{k,i})_{uu}. \end{aligned} \tag{120}\]

These are derivatives of the discrete transition, not unscaled continuous acceleration derivatives. Tassa, Mansard and Todorov’s DDP treatment describes this distinction between full DDP and its Gauss–Newton variant, as well as the need to treat control limits within the local optimization.

17.3 Solving the Unconstrained Quadratic Subproblem

If $Q_{uu}\succ0$ and the input update is unconstrained,

\[ q_u+Q_{xu}^\top\delta x+Q_{uu}\delta u=0, \tag{121}\]

\[ \delta u=K\delta x+d, \tag{122}\]

\[ K=-Q_{uu}^{-1}Q_{xu}^\top,\qquad d=-Q_{uu}^{-1}q_u. \tag{123}\]

This sign convention also applies to the iLQR update:

\[ \delta u_k=K_k\delta x_k+d_k. \tag{124}\]

Use linear solves or factorizations in code, not an explicit matrix inverse. Substitution yields

\[ \begin{aligned} V_k^{\rm local}(\delta x) &=q^0-\tfrac12q_u^\top Q_{uu}^{-1}q_u\\ &\quad+v_k^\top\delta x+\tfrac12\delta x^\top P_k\delta x, \end{aligned} \tag{125}\]

\[ \begin{aligned} v_k&=q_x-Q_{xu}Q_{uu}^{-1}q_u,\\ P_k&=Q_{xx}-Q_{xu}Q_{uu}^{-1}Q_{xu}^\top. \end{aligned} \tag{126}\]

The local quadratic problem is solved exactly in exact arithmetic under its hypotheses. That does not make its approximation to the nonlinear problem exact or establish a global optimum.

17.4 The iLQR Approximation

Dropping only the value-gradient-weighted dynamics Hessians gives

\[ \begin{aligned} Q_{xx}&=L_{xx}+A^\top P^+A,\\ Q_{xu}&=L_{xu}+A^\top P^+B,\\ Q_{uu}&=L_{uu}+B^\top P^+B. \end{aligned} \tag{127}\]

A scalar check makes the omitted terms visible. Let $F(x,u)=x^2+xu+u^2$, $L(x,u)=\tfrac12x^2+\tfrac12u^2+0.2xu$ and $V^+(y)=\tfrac12y^2$. At $(x,u)=(0.3,-0.2)$, $F=0.07$, $DF=(0.4,-0.1)$ and $D^2F=\begin{pmatrix}2&1\\1&2\end{pmatrix}$. The exact action-value Hessian is $\begin{pmatrix}1.30&0.23\\0.23&1.15\end{pmatrix}$. iLQR omits $0.07D^2F$, while retaining the state and cross curvature of $L$. Independent finite differences verify the full Hessian.

17.5 Complete Iteration and Verification

Initialize the nominal states by rolling out chosen inputs from the fixed initial state; initialize $v_N=\phi_x$ and $P_N=\phi_{xx}$. During the backward pass, compute the chosen derivative model and check the input Hessian before factorization. If regularization is needed, state which Hessian or value model is changed and recompute consistently.

During the forward pass, use $u_k^{\rm new}=\bar u_k+\alpha d_k+K_k(x_k^{\rm new}-\bar x_k)$ and roll out the nonlinear transition from the original initial state. A line search compares actual cost reduction with the local prediction. Scaling the feedforward step leaves the feedback correction responding to the newly generated states.

With active input limits, solve the appropriate constrained quadratic subproblem; arbitrary clipping can fail when the input Hessian couples components. A rollout satisfies its discrete dynamics, but it need not satisfy state constraints, continuous-time limits between nodes or the physical plant. Report stationarity, feasibility, accepted step size and cost reduction. A tiny improvement or an iteration limit can indicate stagnation.

17.6 Convergence Requires More Than Positive Input Curvature

Positive $Q_{uu}$ makes a particular local control minimization well posed. It does not alone prove convergence of the nonlinear iteration. Local quadratic Newton convergence requires a nonsingular reduced Hessian, sufficient derivative regularity, a close initialization, and eventually accepted full steps. DDP convergence results depend on the particular algorithm and assumptions. iLQR’s Hessian approximation does not generally give superlinear convergence. See the worked control chapter for the precise distinction.

Use nonlinear rollout, a sufficient-decrease acceptance rule, bounded line search, regularization, and a declared failure exit. If regularization changes the computed policy, evaluate that actual policy in the quadratic model rather than substituting an unregularized Schur-complement formula indiscriminately. Constraints require a compatible constrained subproblem.

No benchmark package establishes a universal 5–20 iteration count here. Dense Riccati algebra includes state-control cross products as well as matrix factorizations; derivative and rollout costs depend on the model and implementation. Report timings together with dimensions, horizon, tolerances, derivative method, and achieved residuals.

17.7 Constrained MPC Extension

17.7.1 Handling Control Bounds

For $u_{\min} \leq u \leq u_{\max}$:

Constrained local solve:

Minimize the local quadratic model subject to bounds on the actual updated input. With a positive-definite input Hessian, this is a convex box-constrained quadratic program. An active-set method must account for coupling between free and saturated components. Componentwise clipping is a solution only in special separable cases; the coupled example above shows why.

Quadratic Penalty:

An input-violation penalty can use $\rho \|\max(0, u-u_{\max})\|^2+\rho \|\max(0,u_{\min}-u)\|^2$. This is not an augmented Lagrangian, which also uses multiplier terms and updates. A finite penalty need not enforce exact feasibility, and increasing $\rho$ can worsen conditioning. Check constraint residuals explicitly. Likewise, clipping an unconstrained update is not equivalent to solving the box-constrained quadratic subproblem.

17.7.2 State Constraints

For path constraints $g(x_k) \leq 0$:

Barrier methods:

\[ L_{\text{aug}}(x, u) = L(x, u) - \mu \sum_i \log(-g_i(x)) \tag{128}\]

A logarithmic barrier requires a strictly feasible interior. Convergence of a sequence as $\mu\to0$ requires appropriately solved subproblems and regularity assumptions; it is not a general global-optimality guarantee for a nonlinear problem.

Sequential quadratic programming (SQP): Linearize constraints, solve QP subproblem per iteration.

17.8 Model Predictive Control (MPC)

17.8.1 Receding Horizon Strategy

At each time $t$:

Measure current state $x(t)$
Solve finite-horizon optimal control (via DDP/iLQR): \[\min_{u(\cdot)} \int_t^{t+T} L(x, u) \, d\tau + \phi(x(t+T)) \tag{129}\]
Apply first control $u^*(t)$ (open-loop for $\Delta t$)
Shift horizon to $[t+\Delta t, t+T+\Delta t]$, repeat

Receding-horizon replanning can respond to the following effects, subject to estimation quality, feasibility, and robustness of the chosen design:

Model mismatch
Disturbances
Unmodeled dynamics

Replanning uses updated state estimates and a chosen predictive model. Its benefit must be assessed with sensing delay, optimization time, admissibility and model uncertainty included.

17.8.2 Stability Guarantees

MPC Stability Theorem

For the sampled model $x^+=F(x,u)$, specify a terminal set $X_f$ and terminal controller $k_f$. Require that $k_f$ respects the constraints, maps $X_f$ into itself, and satisfies

\[ \phi(F(x,k_f(x)))-\phi(x)\leq-\ell_d(x,k_f(x)),\qquad x\in X_f, \]

where $\ell_d$ is the stage cost over a sampling interval. With initial feasibility, exact nominal dynamics, and an attained finite-horizon minimum at each step, shifting the previous optimal sequence and appending $k_f$ supplies a feasible next candidate and the value decrease $V_N(x^+)-V_N(x)\leq-\ell_d(x,u^*)$. A suboptimal solver needs a compatible certified-decrease argument; a small optimization residual alone does not prove this inequality.

Positive-definite state-cost lower bounds, appropriate local upper bounds and continuity of $V_N$, and admissible bounded trajectories then yield the corresponding stability and attraction result on the certified feasible domain. A nonnegative terminal cost alone is not a Lyapunov certificate, and merely choosing a long horizon does not replace these conditions. Disturbances and model mismatch need a separate robust or practical-stability argument.

17.9 Computational Complexity and Timing Budgets

Dense Riccati algebra involves terms of order $N(n^3+n^2m+nm^2+m^3)$, with opportunities to exploit structure. This excludes model-specific derivative evaluation, full-DDP dynamics Hessians, factorization retries, and constrained subproblems. A forward rollout costs $N$ transition evaluations per line-search trial; a general nonlinear transition does not universally cost $O(n)$. Finite differences and automatic differentiation also have different costs depending on the computational graph.

At a control rate of 100 Hz, the available sampling interval is 10 ms. That arithmetic is a budget, not evidence that a specified nonlinear optimization meets it. Record hardware, horizon, dimensions, derivative method, worst-case latency, feasibility, and achieved optimization residuals. Warm starts and iteration limits can reduce runtime, but a time limit requires an admissible fallback policy and does not certify convergence. No universal iteration count or real-time benchmark is established here.

17.10 Robust and Stochastic Extensions

17.10.1 Robust MPC

For uncertain dynamics $\dot{x} = f(x, u, w)$ with bounded $\|w\| \leq \epsilon$:

Min-max formulation:

\[ \min_{u(\cdot)} \max_{\|w\| \leq \epsilon} J(x, u, w) \tag{130}\]

Tube MPC: Specify a nominal trajectory and a feedback policy whose error set remains bounded under every modeled disturbance; tighten constraints accordingly. An open-loop input sequence and a causal feedback policy define different min-max problems. The uncertainty set, information pattern and invariant-error argument are part of the formulation.

17.10.2 Stochastic MPC

For an Itô model $dx=f(x,u)dt+\sigma(x,u)dW_t$ with Brownian motion $W_t$, admissible controls must be adapted to the available information. One possible risk-sensitive objective is:

Risk-sensitive cost:

\[ \min_{u(\cdot)} \mathbb{E}\left[ \int_0^T L(x, u) \, dt + \phi(x(T)) \right] + \lambda \text{Var}[J] \tag{131}\]

Certainty equivalence has precise validity in standard linear-quadratic-Gaussian settings with their separation assumptions. It is not a general robustness theorem for a nonlinear plant. In general, $\mathbb E[f(x,u)]\ne f(\mathbb E[x],\mathbb E[u])$, and even small noise can matter near constraints, unstable modes, or mode transitions. Evaluate the stochastic or robust formulation actually being used.

17.11 Summary: DDP/iLQR/MPC

Key insights:

DDP/iLQR solve local quadratic subproblems whose relation to the nonlinear objective depends on the retained derivatives
Riccati recursion propagates cost-to-go backward through tangent spaces
Feedback gains $K_k$ act on local state deviations within the chosen quadratic model
MPC replans; robustness requires additional design assumptions and verification
Real-time feasibility depends on measured computation, problem structure, and an admissible fallback

These methods use linear algebra and local derivatives in different ways. Their convergence, stability, and performance require the separate assumptions stated above.

18 Case Studies and Worked Examples

These examples test different parts of the argument: exact optimization of a stated linear model, configuration-dependent mechanical coupling, and instantaneous underactuation with controllability over time. Reproducible mathematical results are identified explicitly. None is an experimentally validated golf, spacecraft or flight controller.

18.1 Case Study 1: Spacecraft Rendezvous

18.1.1 Model and Coordinate Convention

Use radial $r_x$, along-track $r_y$ and orbit-normal $r_z$ coordinates in the rotating frame of a circular reference orbit. For small relative separation compared with orbital radius, the Clohessy–Wiltshire approximation is

\[ \begin{aligned} \ddot r_x-2n\dot r_y-3n^2r_x&=u_x/m,\\ \ddot r_y+2n\dot r_x&=u_y/m,\\ \ddot r_z+n^2r_z&=u_z/m. \end{aligned} \tag{132}\]

Here $n$ is the circular reference’s orbital rate and $u$ is force in newtons. MIT’s relative-motion derivation obtains the linear equations by neglecting higher-order relative-position terms. A constant-coefficient linear model is exact as a mathematical model; its approximation to nonlinear orbital mechanics remains limited by these assumptions.

Let $x=(r,\dot r)$, $\dot x=Ax+Bu$, with constant mass. Choose a soft terminal penalty and squared-force cost:

\[ J=\tfrac12x(T)^\top Q_fx(T) +\tfrac12\int_0^T u^\top Ru\,dt. \tag{133}\]

This objective is not fuel consumption and does not impose exact rendezvous. Weights must include appropriate units or arise from declared nondimensional scales. For the numerical example, evaluate positions in metres, velocities in metres per second, forces in newtons and time in seconds; define the numerical objective with $Q_f=\operatorname{diag}(10^4,10^4,10^4,10^2,10^2,10^2)$ and $R=I$. Thus these coefficients weight numerical quantities in those specified units; changing units requires transforming the weights.

18.1.2 Reproducible Sampled LQR Result

Set $n=0.001$ rad/s, $m=100$ kg, $T=1000$ s and $N=50$ held-input intervals of $h=20$ s. Obtain the exact sampled matrices from

\[ \exp\!\left[h\begin{pmatrix}A&B\\0&0\end{pmatrix}\right] =\begin{pmatrix}F&G\\0&I\end{pmatrix}. \]

The discrete running weight is $R_d=hR$. Omitting $h$ changes the objective relative to the terminal penalty. With no running state penalty, the backward recursion is

\[ \begin{aligned} P_N&=Q_f,\\ S_k&=R_d+G^\top P_{k+1}G,\\ K_k&=-S_k^{-1}G^\top P_{k+1}F,\\ P_k&=F^\top P_{k+1}F+F^\top P_{k+1}GK_k. \end{aligned} \]

Roll out $u_k=K_kx_k$, $x_{k+1}=Fx_k+Gu_k$ from the actual initial state. This is one Riccati backward pass and one forward rollout, not an iterative nonlinear solver. The following complete calculation uses NumPy and SciPy:

import numpy as np
from scipy.linalg import expm


def sampled_cw_example() -> tuple[np.ndarray, np.ndarray, float]:
    """Solve the stated unconstrained, held-input quadratic problem."""
    rate, mass, horizon, intervals = 0.001, 100.0, 1000.0, 50
    step = horizon / intervals
    a = np.zeros((6, 6))
    a[:3, 3:] = np.eye(3)
    a[3, 0], a[3, 4] = 3 * rate**2, 2 * rate
    a[4, 3], a[5, 2] = -2 * rate, -rate**2
    b = np.zeros((6, 3))
    b[3:] = np.eye(3) / mass
    generator = np.zeros((9, 9))
    generator[:6, :6], generator[:6, 6:] = a, b
    transition = expm(step * generator)
    f, g = transition[:6, :6], transition[:6, 6:]
    terminal = np.diag([1e4, 1e4, 1e4, 1e2, 1e2, 1e2])
    running = step * np.eye(3)
    p = terminal.copy()
    gains = []
    for _ in range(intervals):
        gain = -np.linalg.solve(running + g.T @ p @ g, g.T @ p @ f)
        gains.append(gain)
        p = f.T @ p @ f + f.T @ p @ g @ gain
        p = (p + p.T) / 2
    states = [np.array([100.0, 50.0, 20.0, 0.01, 0.02, 0.0])]
    controls = []
    for gain in reversed(gains):
        control = gain @ states[-1]
        controls.append(control)
        states.append(f @ states[-1] + g @ control)
    x, u = np.asarray(states), np.asarray(controls)
    cost = 0.5 * x[-1] @ terminal @ x[-1]
    cost += 0.5 * sum(control @ running @ control for control in u)
    return x, u, float(cost)

The reproduced terminal position is approximately $(1.66617\times10^{-6},3.04859\times10^{-7},1.48660\times10^{-7})$ m. Its terminal velocity is approximately $(-0.0507127,0.0117398,-0.00839610)$ m/s. The objective is 1.22796775, peak force magnitude 0.09360394 N and accumulated force magnitude 39.39750874 N·s. An independent condensed quadratic-program calculation agrees with the Riccati control sequence to approximately $10^{-8}$ N; exact sampled propagation also agrees with independently integrated continuous dynamics.

The very small position error coexists with substantial terminal velocity. That is the actual tradeoff of these soft weights, not successful docking. Exact endpoint constraints, thrust limits, collision avoidance, actuator allocation and nonlinear orbital-model verification would define a different problem.

For an ideal thruster at fixed specific impulse $I_{\rm sp}$, propellant relates to delivered impulse through $m_p=I_{\rm total}/(g_0I_{\rm sp})$, as follows from NASA’s specific-impulse definition. For example, 12.4 N·s at 300 s corresponds to 0.00421483 kg, using $g_0=9.80665$ m/s². This separate unit check is not the impulse of the simulation above. Real thruster allocation can make the sum of individual thrust magnitudes exceed the norm of the resultant, and squared force is a different integral again.

18.2 Case Study 2: Planar Robot Arm

18.2.1 Placement Determines Both Gravity and Kinematics

Choose two massless links of length 1 m, with a 1 kg point mass at each distal end. This is not the inertia model of two uniform massive rods. Let $q_1$ be measured counterclockwise from the horizontal world $+x$ axis and $q_2$ be the relative elbow angle. Then

\[ \begin{aligned} p_1&=(\cos q_1,\sin q_1),\\ p_2&=p_1+(\cos(q_1+q_2),\sin(q_1+q_2)). \end{aligned} \]

Writing

\[ M(q)\ddot q+C(q,\dot q)\dot q+G(q)=\tau, \tag{134}\]

the chosen masses and lengths give

\[ M(q)= \begin{pmatrix} 3+2\cos q_2&1+\cos q_2\\ 1+\cos q_2&1 \end{pmatrix}, \tag{135}\]

in kg·m². With $a=-\sin q_2$ in the corresponding mass-length units,

\[ C(q,\dot q)= \begin{pmatrix} a\dot q_2&a(\dot q_1+\dot q_2)\\ -a\dot q_1&0 \end{pmatrix}. \tag{136}\]

The potential is $U=g[2\sin q_1+\sin(q_1+q_2)]$ in joules with these SI parameter values, so

\[ G(q)=g \begin{pmatrix} 2\cos q_1+\cos(q_1+q_2)\\ \cos(q_1+q_2) \end{pmatrix}. \tag{137}\]

Gravity contains cosines for angles measured from the horizontal. A sine formula would require a different angular origin and matching kinematics.

Independent checks construct $M$ from the two point-mass Jacobians, differentiate $U$ numerically and verify $\dot q^\top(\dot M-2C)\dot q=0$. Thus the mechanical energy satisfies $\dot E=\dot q^\top\tau$. These checks link placement, inertia, gravity and power; checking only that $M$ is positive definite would miss the inconsistent gravity convention.

18.2.2 Inverse Kinematics Before Trajectory Optimization

For a target $(x,y)$ with general positive link lengths, define $c_2=(x^2+y^2-L_1^2-L_2^2)/(2L_1L_2)$. Reachability requires $|L_1-L_2|\leq\sqrt{x^2+y^2}\leq L_1+L_2$. Reject targets outside this annulus apart from a stated numerical tolerance; clipping an arbitrary unreachable target into $[-1,1]$ invents a different target.

For the branch $q_2=+\arccos c_2$, take $q_1=\operatorname{atan2}(y,x)-\operatorname{atan2}(L_2\sin q_2,L_1+L_2\cos q_2)$. The negative branch gives the other elbow configuration when distinct. At a singular folded configuration, additional care is needed because the shoulder angle can be indeterminate.

With unit lengths, the checked angles are:

Target $(0.5,0.5)$ m: $q\approx(-24.29518895^\circ,138.59037789^\circ)$.
Target $(0.3,0.8)$ m: $q\approx(4.73396255^\circ,129.41998447^\circ)$.

By contrast, $(30^\circ,60^\circ)$ places the end at $(\sqrt3/2,1.5)$ m. Forward substitution verifies each result. These are kinematic endpoints, not a computed dynamically feasible trajectory or a convergence benchmark.

18.2.3 What Coupling Establishes

The off-diagonal entries of $M$ mean that a torque at one joint can accelerate both joints. Their dependence on posture changes that instantaneous mapping. It does not mean the tangent-space dimension changes, and matrix entries do not define a universal ordering of configurations from least to greatest inertia.

DDP can use derivatives along a changing nominal motion, while a fixed linear controller can still work within a verified operating region. Neither possibility is decided solely by observing that $M(q)$ varies. For a golf model, anatomical constraints, shaft modes, ground reactions and available actuator directions add further structure; the fully actuated two-link arm does not establish a universal optimal sequencing rule.

18.3 Case Study 3: Quadrotor Stabilization

18.3.1 Hover in a Consistent Frame

Let world and body $e_3$ point upward at level attitude, and let $R$ map body vectors to world vectors. Rotor thrusts $u_i\geq0$ act along body $+e_3$. Let $q$ be a unit, scalar-first quaternion for $R$ and $\omega$ be body angular velocity. Then

\[ \begin{aligned} \dot p&=v,\\ \dot v&=-g e_3+\frac{\sum_i u_i}{m}R e_3,\\ \dot q&=\tfrac12q\otimes(0,\omega),\\ \dot\omega&=J^{-1}(\tau-\omega\times J\omega). \end{aligned} \tag{138}\]

The pure quaternion $(0,\omega)$ is essential to the quaternion product. This continuous equation preserves $\|q\|=1$; a numerical method must also control its constraint error.

Place rotors at body locations $(L,0,0)$, $(0,L,0)$, $(-L,0,0)$ and $(0,-L,0)$. A consistent alternating-spin allocation is

\[ \begin{aligned} \tau_x&=L(u_2-u_4),\\ \tau_y&=L(u_3-u_1),\\ \tau_z&=c(u_1-u_2+u_3-u_4). \end{aligned} \]

Here $u_i$ is thrust in newtons, $L$ is in metres and $c$ is a torque-to-thrust ratio in metres. These inputs are not raw motor speed commands: if $u_i=k_f\Omega_i^2$, the model is affine in thrust but nonlinear in $\Omega_i$.

At $R=I$, zero rates and any fixed position, $u_i=mg/4$ balances gravity and produces zero moment. Adding gravity and upward thrust with the same sign would not admit this hover equilibrium.

18.3.2 Twelve Local Errors for a Thirteen-Component State

The embedded vector $(p,v,q,\omega)$ has 13 entries but lies on a 12-dimensional state manifold. Near a fixed hover reference, use the stacked error

\[ e=\begin{pmatrix} p-p_{\rm ref}\\ v\\ \operatorname{Log}(R_{\rm ref}^\top R)^\vee\\ \omega \end{pmatrix}. \]

The logarithm is a local rotation coordinate, with branch restrictions away from the identity; it is not a global nonsingular three-parameter attitude representation.

A $12\times12$ quadratic weight acts on this 12-component error, not on the raw 13-component subtraction. The construction also respects that $q$ and $-q$ represent the same attitude. A running cost can penalize $e$ and $\delta u=u-u_{\rm hover}$ so the desired steady operation has zero deviation cost. Moving references require consistent angular-rate transport and tracking errors.

18.3.3 Instantaneous Rank and Controllability Over Time

For $m=1.2$ kg, $L=0.2$ m, $c=0.01$ m and $J=\operatorname{diag}(0.02,0.02,0.04)$ kg·m², the linearized input matrix has rank four. At level hover, $\delta\dot v_x=g\,\delta\theta_y$ and $\delta\dot v_y=-g\,\delta\theta_x$: rotor-induced moments first change attitude, and attitude then changes horizontal acceleration.

The controllability matrix $[B,AB,\ldots,A^{11}B]$ has rank 12 for this example. This does not create twelve independent instantaneous inputs. It shows that the linearized state can be controlled over time near an interior hover input, subject to the local-model assumptions. Positive thrust bounds and available thrust margins still constrain finite maneuvers.

At 40 Hz, 50 intervals span 1.25 s; a two-second horizon requires 80 intervals. These are sampling facts. Settling time, centimetre tracking accuracy and robustness percentages require an implemented controller, declared perturbations and actual verification; no such performance claim is established here.

18.4 What the Examples Contribute to Swing Analysis

The orbital example shows why an exactly solved quadratic objective can still produce an unacceptable terminal event. The arm shows that kinematics, inertia, gravity and energy must share one physical placement. The quadrotor shows why instantaneous input rank and the ability to influence a future task are different questions.

Together they motivate a disciplined swing argument: derive a consistent mechanical model, specify the available inputs and intervention, propagate the resulting motion, and evaluate the actual contact objective with its constraints and uncertainty. Local derivatives make parts of this process tractable; they do not replace the physical and experimental checks.

18.5 Reproduction Scope

The companion control and case-study checks reproduce the Hamiltonian signs, constrained quadratic example, integration orders, DDP curvature, sampled orbital calculation, arm mechanics and hover rank calculations. They are verification examples rather than a released nonlinear-control package or empirical validation.

19 Conclusion and Future Directions

19.1 Three Levels of Mechanical Reasoning

The useful result is a hierarchy of statements with different scopes.

An immediate acceleration change. In a model affine in the declared force input, the increment is $M^{-1}\mathcal B\,\Delta u$. State, time, model, parameters and constraint mode remain fixed.
A small change in a later output. The derivative of the solution/output map propagates through the variational equation. Reference motion, perturbation definition, policy and endpoint convention remain fixed.
A finite intervention. A new nonlinear trajectory determines its own task and event outputs. The intervention and comparison protocol are specified, while the state evolves independently.

A tangent space supplies a place for velocity vectors. Differentiation supplies a linear map between tangent spaces. Neither fact makes the underlying vector field or its finite flow linear. First variations superpose because they obey one linear equation along a common reference, not because all nonlinear motion is a sequence of exactly linear physical systems.

19.2 How the Swing’s Components Interact

At one posture, an applied joint torque can accelerate several coordinates through the inverse inertia matrix. Shaft deflection changes the geometry and energy stores entering that matrix and the other loads. Subsequent posture and velocity changes alter both the input gain and the reference transport. Finally, clubhead kinematics and contact timing select which of those state changes matter for the strike.

This sequence connects local mechanical coupling to an eventual task outcome without identifying them. A large instantaneous acceleration contribution can have a small endpoint effect. Positive input work can accompany a decrease in a particular speed. A small interaction residual between two perturbations can coexist with substantial model bias or a large single-trajectory approximation error.

Input definition also matters. Replaying a torque history, preserving a feedback policy, changing stiffness, changing excitation in an activation model, and reoptimizing after an equipment change are different experiments. A result from one should not be described as the result of another.

19.3 What Linear Control Can Establish

For an equilibrium with the required smoothness, a Hurwitz closed-loop Jacobian can support a local nonlinear stability theorem. LQR additionally solves a specified linear-quadratic optimization problem under its stabilizability and detectability conditions. Finite-horizon trajectory control, DDP/iLQR and MPC need their own conditions on feasible references, objectives, constraints, numerical rollout and closed-loop execution.

An exact derivative is valuable input to those methods. It does not by itself establish global optimality, a region of attraction, a guaranteed contact result or a biological explanation. A certificate must state the model, domain, disturbance assumptions and physical output it controls.

19.4 A Verification Protocol for the Argument

For a proposed golf application, first specify the state, coordinate convention, actuator inputs, constraints and measured outputs. Demonstrate dimensional consistency and force/energy balances. Then test derivative code against independent finite perturbations on smooth intervals. Refine the integration tolerance and perturbation amplitude separately.

Compare outcomes at the same declared time or event. Include event-time and reset sensitivity when appropriate. Test the actual proposed feedback or input replay, retain correlated uncertainty, and report the output scale and tolerance used to judge a discrepancy. Verify a finite nonlinear intervention directly before presenting its first-order estimate as quantitatively adequate.

Finally, distinguish numerical verification from empirical validation. An internally consistent example can test an equation or algorithm. Evidence for a golfer requires identified parameters, measured loads and motion, uncertainty analysis, and prediction on data not used to fit the model. Local mathematical correctness is necessary for that investigation, but cannot substitute for it.

19.5 Open Technical Questions

The next research steps are concrete: identify which input directions can be estimated from available measurements; determine the domain over which finite predictions meet task tolerances; and test how contact transitions, flexible modes, actuator limits and feedback alter those conclusions. A proposed geometric residual bound must distinguish coordinate Hessians, dynamic transport and intrinsic curvature. A stochastic extension needs a specified stochastic differential equation and calculus convention; Brownian paths do not have ordinary instantaneous velocities to which the deterministic argument can simply be applied.

These questions do not have universal answers based on swing duration or the dimension of a tangent space. They can be investigated through explicit models, controlled comparisons and evidence.

In Plain Language

A snapshot can tell us how an extra torque changes acceleration at the current posture. Predicting the rest of the swing requires following how that change alters the body’s motion, the shaft and the later effect of every load. A small-change calculation is a precise way to study that connection near a reference swing. A larger change needs a new simulation and, for claims about real golfers, measurements that test its prediction.

References

Abraham, R., & Marsden, J. E. (1978). Foundations of Mechanics (2nd ed.). Benjamin/Cummings.
Arnold, V. I. (1989). Mathematical Methods of Classical Mechanics (2nd ed.). Springer.
Bryson, A. E., & Ho, Y.-C. (1975). Applied Optimal Control. Hemisphere Publishing.
Jacobson, D. H., & Mayne, D. Q. (1970). Differential Dynamic Programming. Elsevier.
Khalil, H. K. (2002). Nonlinear Systems (3rd ed.). Prentice Hall.
Lee, J. M. (2012). Introduction to Smooth Manifolds (2nd ed.). Springer.
Liberzon, D. (2003). Switching in Systems and Control. Birkhäuser.
Mayne, D. Q., et al. (2000). Constrained model predictive control: Stability and optimality. Automatica, 36(6), 789-814.
Sastry, S. (1999). Nonlinear Systems: Analysis, Stability, and Control. Springer.
Tassa, Y., et al. (2012). Synthesis and stabilization of complex behaviors through online trajectory optimization. In Proc. IEEE/RSJ IROS (pp. 4906-4913).

Appendices

20 Appendix A: Mathematical Prerequisites

A smooth manifold is described locally by coordinates with smooth invertible transition maps. $T_x\mathcal M$ contains curve velocities at one point; $T\mathcal M$ collects these spaces over the manifold. A covector is a linear functional on a tangent space. Generalized force and the differential of a scalar cost are covectors, whereas velocity and a first state variation are vectors.

A smooth map has a derivative between the appropriate tangent spaces. Coordinates represent this derivative by a Jacobian. A metric supplies lengths and identifies vectors with covectors, while a connection supplies a rule for differentiating vector fields along curves. Neither structure is required merely to define a tangent space or the derivative of a smooth map. If used, each must be specified.

For rotations, use a local chart or a declared Lie-group representation with its frame convention. A coordinate singularity is not automatically a physical instability. For a mechanical state, distinguish configuration from its velocity state, and distinguish both from a clubhead output. These definitions are developed with physical examples in the configuration-manifold chapter.

21 Appendix B: Proofs and Derivations

This appendix supplies explicit bounds and convergence statements used in the main text.

22 B.1 A Finite-Difference Bound From the Endpoint Map

Let $F(a,b)=x(t_1;\bar u+a v_1+b v_2)$ use the same initial state for every pair of amplitudes. Assume a $C^2$ solution map exists on a rectangle containing $0\leq a,b\leq\epsilon$. The fundamental theorem of calculus, applied twice componentwise, gives the exact identity

\[ r=F(\epsilon,\epsilon)-F(\epsilon,0)-F(0,\epsilon)+F(0,0) =\int_0^\epsilon\int_0^\epsilon F_{ab}(a,b)\,db\,da. \]

Consequently,

\[ \|r\|\leq\epsilon^2\sup_{a,b\in[0,\epsilon]}\|F_{ab}(a,b)\|. \]

The mixed second-variation equation in the main text computes $F_{ab}$, using derivatives along each corresponding trajectory. A bound on the state Hessian alone omits mixed state-input terms and nonlinear input effects. For control-affine dynamics $f_0(x)+G(x)u$, the term $D G(x)[\xi]v$ remains even though the pure input Hessian vanishes. Existence of one nominal trajectory is also insufficient: all perturbations used by the bound must remain in the smooth model domain.

23 B.2 Intrinsic Curvature Is a Different Quantity

For a connection $\nabla$, the curvature identity is

\[ R(v,w)\xi=\nabla_v\nabla_w\xi-\nabla_w\nabla_v\xi-\nabla_{[v,w]}\xi. \]

The final Lie-bracket term vanishes for commuting coordinate vector fields, but cannot generally be omitted. Curvature can govern parallel-transport holonomy. That geometric operation is not the same as comparing four forced trajectories of an arbitrary vector field.

The counterexample $\dot x=u^2$ lives on the flat Euclidean line, where $R=0$, yet its input interaction is nonzero. It disproves a general bound of this interaction solely by sectional curvature. Coordinate-dependent field derivatives, a specified connection, and intrinsic manifold curvature must be distinguished before proposing a geometric interpretation.

24 B.3 Integration Error and Linearization Error

For explicit Euler applied to a sufficiently smooth initial-value problem over a fixed finite interval, bounded local truncation error $O(h^2)$ and a suitable Lipschitz stability bound yield global integration error $O(h)$. The constant depends on the interval and growth of perturbations. Lipschitz continuity alone does not supply the higher differentiability needed for every claimed high-order method.

If the integrated model is the first-variation approximation rather than the full nonlinear system, reducing $h$ converges to that approximate model. It does not remove finite-perturbation error. For a single smooth input direction, write $x_\epsilon=\bar x+\epsilon\xi+\tfrac12\epsilon^2\eta+o(\epsilon^2)$. A computed first-variation endpoint with a uniform Euler error bound satisfies

\[ \|x_\epsilon(t_1)-x_{\mathrm{linear},h}(t_1)\| \leq C_{\mathrm{disc}}h+\tfrac12\epsilon^2\|\eta(t_1)\|+o(\epsilon^2). \]

The self-variation $\eta$ includes state, input, and mixed derivatives and their transport. This distinguishes convergence of the numerical solver from fidelity of the linearized physical prediction. Assess the two with separate step-size and perturbation-amplitude refinements.

25 Appendix C: Notation Guide and Glossary

26 Primary Notation

Symbol	Meaning
$x \in \mathbb{R}^n$	State vector
$u \in \mathbb{R}^m$	Control input
$f: \mathbb{R}^n \times \mathbb{R}^m \to \mathbb{R}^n$	Vector field (dynamics)
$A = \frac{\partial f}{\partial x}$	State Jacobian
$B = \frac{\partial f}{\partial u}$	Input Jacobian
$\delta x$, $\delta u$	Infinitesimal perturbations
$\Phi(t_1, t_0)$	State transition operator
$\lambda(t)$	Costate (adjoint variable)
$H(x,u,\lambda,t)$	Hamiltonian
$L(x,u,t)$	Optimization running cost; distinct from a mechanical Lagrangian
$r$	Residual (failure of superposition)

27 Glossary

Costate ($\lambda$): A covector propagated by an adjoint equation for a specified cost and comparison protocol. An open-loop future-cost sensitivity need not equal the gradient of an optimized value function.

Fréchet Derivative: A bounded linear map whose remainder is little-o of the perturbation norm. Existence of derivatives along separate lines is insufficient.

Jacobian: Coordinate matrix of a differentiable map’s derivative. A state-rate Jacobian along a moving coordinate transformation includes the associated basis-change term.

Manifold: Space that locally resembles Euclidean space; may be curved globally.

Residual: Deviation from additivity when applying multiple perturbations; $O(\epsilon^2)$ interaction under stated smoothness and fixed-horizon assumptions.

State Transition Operator ($\Phi$): Linear map transporting perturbations forward in time via variational dynamics.

Tangent Space ($T_x \mathcal{M}$): Linear space of velocity vectors at point $x$ on manifold $\mathcal{M}$.

Variational Equation: The exact evolution equation for the derivative of a trajectory family with respect to perturbation amplitude, along one common reference. Its finite-amplitude use is an approximation.

28 Appendix D: Code Repository Index

The companion regression file test_tangent_foundation_rigor.py independently checks the differentiability counterexamples, coordinate transformation, pendulum energy, ordered transport, pulse kernel, functional/adjoint sensitivity and combined-motion work example. The existing residual checks check the evolved baseline and mixed-variation transport. These are verification examples, not a released general-purpose nonlinear-control package. The control and case-study checks independently verify Part III’s optimization, discretization, orbital, arm and hover examples. The sampled orbital calculation is reproducible; no empirical flight or golf performance is claimed.

End of Thesis

Cite this page

Citation

BibTeX citation: Download .bib

@online{olson2026,
  author = {Olson, Dieter},
  title = {Full {Reference} {Manuscript:} {Exact} {Instantaneous}
    {Superposition} in {Nonlinear} {Dynamics}},
  date = {2026-01-01},
  url = {https://affinedrift.com/articles/tangent-hyperplane-articles/Tangent_Hyperplanes_Unified_Thesis.html},
  langid = {en}
}

For attribution, please cite this work as:

Olson, Dieter. 2026. “Full Reference Manuscript: Exact Instantaneous Superposition in Nonlinear Dynamics.” January 1. https://affinedrift.com/articles/tangent-hyperplane-articles/Tangent_Hyperplanes_Unified_Thesis.html.

--- title: "Full Reference Manuscript: Exact Instantaneous Superposition in Nonlinear Dynamics" description: "Reference manuscript for the tangent-space series. Use the compact seven-part series as the canonical reading path." subtitle: "Reference Derivations for Tangent Space Linearization and Integral Accumulation" author: "Dieter Olson" date: "2026-01-01" date-source: "initial-publication-record" format: pdf: toc: true toc-title: "Contents" toc-depth: 3 number-sections: true pdf-engine: xelatex html: toc: true toc-depth: 4 number-sections: true code-fold: true fig-cap-location: bottom execute: echo: false bibliography: references.bib categories: - tangent-space --- ::: {.callout-warning} ### Reference Status This is the single-file full reference manuscript for the tangent-space material. It is retained for derivations, appendices, and historical traceability. The canonical rendered reading path is the compact seven-part [Tangent-Space Series](../tangent-hyperplanes-series/part-1-geometry.html), which now carries the contraction, hybrid, and residual-aware topics in shorter source-backed form. ::: ::: {.callout-note appearance="minimal"} ### Thesis Abstract A golf swing couples changing posture, segment momentum, shaft deformation, applied torque and contact. Several useful linear relations survive inside this nonlinear motion, but they answer different questions. A differentiable model has a linear derivative at a specified state and input. A control-affine model additionally has an exact finite input-to-rate increment relation at a fixed state. The first variation of a trajectory satisfies a linear differential equation along a common reference motion. None of these facts makes independently evolving swings additive. Part I develops these distinctions, their geometric meaning, and the regularity needed for error bounds. Part II derives how an entire input perturbation is transported into an endpoint change and how an adjoint weights it by a chosen outcome. Part III develops optimal-control methods subject to their own assumptions. The purpose is a technically usable account of local modeling, with explicit limits on finite predictions and physical interpretation. ::: ::: {.callout-note} ### A Note on "Exactness" The derivative and first-variation equation are exact mathematical objects. Replacing a finite nonlinear response by its first-order expansion is an approximation. Continuous differentiability gives a little-o remainder; a quadratic bound needs additional regularity. Exact fixed-state force addition is a separate mechanical statement about a specified input map. Keep these meanings distinct throughout the [canonical rendered reading path](../tangent-hyperplanes-series/part-1-geometry.html). ::: ## Part I: Geometric Foundations and Local Dynamics {.unnumbered} ## Introduction and Motivation ### The Central Thesis Suppose a model predicts the body's and club's future state from initial conditions and an applied-torque history. Three experiments ask three different questions: 1. At the same posture, velocity and shaft state, how does an extra torque change acceleration now? 2. Along a specified reference swing, how does a small perturbation change the eventual clubhead state to first order? 3. If that torque history is actually changed by a finite amount, what new motion and impact follow? The first can have an exact affine answer. The second is a derivative of the model's solution map. The third generally requires a new nonlinear integration, including any changed event time. Their relationship explains the usefulness of linear methods without extending their conclusions beyond their hypotheses. ### Central Claim For a differentiable coordinate model, the map from state/input perturbation directions to the **first variation of state rate** is linear. If the model is also affine in the declared input, finite input increments produce exactly linear rate increments while state, time and mode remain fixed. Along a common reference, first variations propagate linearly. These are properties of the model and its derivatives; they do not establish biological intent, individual muscle recruitment or optimal swing technique. ### Scope: No Trajectory Superposition For a homogeneous linear system, linear combinations of solutions are solutions. For a linear system with a fixed nonzero forcing, differences about a common particular solution obey the homogeneous equation. Even in the linear case, adding two complete forced solutions would double the common forcing. For a nonlinear system, the finite input-to-trajectory map generally fails to be additive. This can happen on ordinary Euclidean space with its usual flat geometry. Curvature of the state manifold is neither required nor a universal measure of the error. ### Terminology: "Tangent Hyperplane" as Affine Subspace At fixed state $x$, a control-affine model has rates $f_0(x)+G(x)u$. For unrestricted $u\in\mathbb R^m$, this image is the affine subspace $f_0(x)+\operatorname{im}G(x)$ of $T_x\mathcal M$. It is a hyperplane only when its codimension is one. With bounded admissible inputs $\mathcal U(x)$, the actual instantaneous rate set is $f_0(x)+G(x)\mathcal U(x)$, usually only part of that subspace. This set is not the finite-time reachable set. ### Scope and Prerequisites We use a finite-dimensional state chart and one smooth contact/constraint mode. Assume a unique reference solution exists on the stated finite interval, and all compared solutions remain in the model domain. For the initial derivations, take $f$ continuous in time and continuously differentiable in state and input on a neighborhood of those solutions. Piecewise continuous prescribed inputs can be handled on each smooth time segment. More general measurable inputs require corresponding uniform regularity and integrability hypotheses. State and input components should have declared units or nondimensional scales before norms are compared. A configuration $q$ is not the complete dynamic state: velocity, shaft modes and retained actuator states matter. At impacts, saturation switches, friction transitions or changing constraint rank, ordinary smooth derivatives may cease to describe the experiment. A smooth equality-constrained system can be reduced to a regular chart; an active inequality or contact transition needs additional analysis. ## Tangent Spaces and Exact Differentiation ### Nonlinear Dynamics and Smoothness In local coordinates, let $$ \dot x=f(x,u,t),\qquad x\in\mathbb R^n,\quad u\in\mathbb R^m. $$ {#eq-thesis-dynamics} Fix $(\bar x,\bar u,t)$ and combine the perturbations into $h=(h_x,h_u)$. Fréchet differentiability means there is a linear map $D f$ such that $$ f(\bar x+h_x,\bar u+h_u,t) =f(\bar x,\bar u,t)+D f(\bar x,\bar u,t)[h]+r(h), $$ {#eq-frechet} $$ \lim_{\|h\|\to0}\frac{\|r(h)\|}{\|h\|}=0. $$ {#eq-little-o} The limit controls arbitrary directions approaching zero in the chosen norm. It does not say $r(h)=0$ at a nonzero perturbation. Differentiability at a point does not require a continuous derivative in a whole neighborhood; $C^1$ is a useful stronger hypothesis that gives uniform control on compact neighborhoods and supports the flow results used here. ### Continuous Differentiability Does Not Give a Quadratic Bound The scalar function $F(z)=|z|^{3/2}$ is $C^1$ near zero and $F'(0)=0$. Its remainder is $|h|^{3/2}$. Dividing by $|h|$ gives a quantity tending to zero, as required for differentiability. Dividing by $h^2$ gives $|h|^{-1/2}$, which is unbounded. Thus a blanket $O(\|h\|^2)$ claim under only $C^1$ assumptions is false. A locally Lipschitz derivative, often denoted $C^{1,1}$, is sufficient for a quadratic bound. If $\|Df(z+s h)-Df(z)\|\leq L_D s\|h\|$ along the segment, then $$ r(h)=\int_0^1\bigl(Df(z+s h)-Df(z)\bigr)h\,ds, \qquad \|r(h)\|\leq\tfrac12 L_D\|h\|^2. $$ A continuous bounded second derivative on the relevant neighborhood also supplies this bound. The neighborhood, norm, constant and permissible amplitudes are part of the statement, not optional details. ### Directional Derivatives For a Fréchet-differentiable map, each directional derivative agrees with the same linear map: $$ \left.\frac{d}{d\epsilon}f(z+\epsilon h,t)\right|_{\epsilon=0} =Df(z,t)[h]. $$ {#eq-directional} Existence of derivatives along every line does not by itself establish this property. Define $F(a,b)=a^3/(a^2+b^2)$ away from the origin and $F(0,0)=0$. It is continuous at the origin, and its derivative along direction $(a,b)$ is $a^3/(a^2+b^2)$. The directions $(1,0)$ and $(0,1)$ give 1 and 0, while their sum gives $1/2$. These directional derivatives do not constitute a linear derivative. ### The Jacobian Matrices Under differentiability, the coordinate derivative has blocks $$ A(t)=f_x(\bar x(t),\bar u(t),t),\qquad A_{ij}=\frac{\partial f_i}{\partial x_j}, $$ {#eq-jacobian-A} $$ B(t)=f_u(\bar x(t),\bar u(t),t),\qquad B_{ij}=\frac{\partial f_i}{\partial u_j}. $$ {#eq-jacobian-B} The exact first variation of the rate is $A h_x+B h_u$. Both matrices use the same reference state, input and time. A derivative evaluated on one counterfactual branch cannot silently be substituted for that of another. ### A Tangent Space Is Not a Derivative The tangent space $T_x\mathcal M$ is the vector space of velocities of curves through $x$ on a smooth manifold. It exists independently of any particular dynamic law. A smooth map $F:\mathcal M\to\mathcal N$ has a derivative $T_xF:T_x\mathcal M\to T_{F(x)}\mathcal N$. The space and the map are different objects. A vector field assigns $f(x)\in T_x\mathcal M$. Geometrically it is a section $\mathcal M\to T\mathcal M$; its derivative therefore maps into the tangent space of the **tangent bundle**, not automatically back into $T_x\mathcal M$. Local coordinates or an explicitly chosen connection give the familiar component formulas. A state-rate Jacobian along a moving reference is a valid coordinate representation of the variational dynamics; it is not an intrinsic endomorphism independent of that choice. For a family of solutions $x_\epsilon(t)$ with prescribed $u_\epsilon=\bar u+\epsilon v$ and $x_\epsilon(t_0)=\bar x_0+\epsilon\eta_0$, define $$ \eta(t)=\left.\frac{\partial x_\epsilon(t)}{\partial\epsilon}\right|_0, \qquad \dot\eta=A(t)\eta+B(t)v(t),\quad\eta(t_0)=\eta_0. $$ {#eq-tangent-dynamics} This equation follows by differentiating the dynamic equation with respect to the family parameter. $\eta$ is a derivative with respect to amplitude, not the full finite displacement $x_\epsilon-\bar x$. Under the stated regularity on a fixed common interval, the latter is $\epsilon\eta+o(\epsilon)$. ### Coordinate Changes Along a Moving Motion For a smooth invertible coordinate change $y=\psi(x)$, write $T(t)=D\psi(\bar x(t))$. Variations satisfy $\eta_y=T\eta_x$. Differentiating this identity gives $$ A_y=(\dot T+T A_x)T^{-1},\qquad B_y=T B_x. $$ The additional $\dot T$ term is essential away from an equilibrium. For $\dot x=1$ on $x>0$ and $y=x^2$, $A_x=0$, while $\dot y=2\sqrt y$ has $A_y=1/x$. A constant $\eta_x$ corresponds to $\eta_y=2x\eta_x$, which grows as $x$ grows. This is the same physical separation represented in different units and coordinates. At an equilibrium of an autonomous system, $\dot T=0$ and the Jacobians are similar, so their eigenvalues agree. This matters when comparing a joint-space perturbation with a clubhead-space perturbation. A growing coordinate norm need not mean that the physically chosen task error grows at the same rate. ### Example: The Simple Pendulum For a point mass $m$ at length $\ell$, torque $u$, and angle $\theta$ measured from downward vertical, $$ m\ell^2\ddot\theta+mg\ell\sin\theta=u. $$ {#eq-pendulum-second-order} With $x=(\theta,\omega)$, $$ f(x,u)=\begin{pmatrix}\omega\\-(g/\ell)\sin\theta+u/(m\ell^2)\end{pmatrix}. $$ {#eq-pendulum-state-space} The derivative is $$ D f[h_x,h_u]=A h_x+B h_u, $$ {#eq-pendulum-jacobian-def} $$ A=\begin{pmatrix}0&1\\-(g/\ell)\cos\bar\theta&0\end{pmatrix}, \qquad B=\begin{pmatrix}0\\1/(m\ell^2)\end{pmatrix}. $$ {#eq-pendulum-jacobian} Along any reference solution, $$ \dot\eta_\theta=\eta_\omega,\qquad \dot\eta_\omega=-(g/\ell)\cos\bar\theta(t)\eta_\theta +v/(m\ell^2). $$ {#eq-pendulum-linearized} At the unforced downward equilibrium, this becomes $$ \ddot\eta_\theta+(g/\ell)\eta_\theta=0, $$ {#eq-pendulum-harmonic} $$ \lambda_{1,2}=\pm i\sqrt{g/\ell}. $$ {#eq-pendulum-eigenvalues} Pure imaginary eigenvalues do not prove stability of the nonlinear equilibrium by the indirect method. Here a separate energy argument does: $E=\tfrac12m\ell^2\omega^2+mg\ell(1-\cos\theta)$ is conserved and has a strict local minimum at $(0,0)$. Sufficiently small sublevel sets remain in a small neighborhood. The equilibrium is locally Lyapunov stable, but nearby positive-energy motions cannot converge to it without dissipation. For comparison, the planar systems $\dot z=Jz\pm\|z\|^2z$, with $J$ a unit rotation generator, have the same imaginary eigenvalues at the origin. Their radii satisfy $\dot r=\pm r^3$: one decays and the other grows. The derivative alone leaves that distinction unresolved. The pendulum's finite-angle discrepancy can be seen from $$ \sin\theta=\theta-\theta^3/6+O(\theta^5). $$ {#eq-sin-taylor} At this particular equilibrium the first omitted state term is cubic by symmetry. Generic nonlinear models need not have that cancellation. ### Exact Superposition of Input Increments Linearity of $Df$ gives additivity of derivative directions for any differentiable model. Exact finite input additivity needs the stronger structural condition $$ f(x,u,t)=f_0(x,t)+G(x,t)u. $$ At fixed $(x,t)$, with one nominal input $\bar u$, define $\Delta_i=f(x,\bar u+v_i,t)-f(x,\bar u,t)$. Then the combined increment is exactly $\Delta_{12}=G(v_1+v_2)=\Delta_1+\Delta_2$. The baseline rate appears once. This statement can hold for finite increments without making the state dependence linear. For regular mechanical coordinates, $M(q)\ddot q+h(q,\dot q,t)=\mathcal B(q)u$ has this property when $M$ is nonsingular and neither $M$ nor $h$ depends on the declared input. Applied forces are covectors; $M^{-1}\mathcal B$ maps them into acceleration increments. Equality-constraint reactions can preserve affinity after a regular linear solve. A fixed contact mode alone does not guarantee it: a nonlinear actuator law, such as torque proportional to the square of a command, must be inspected in the actual command variable. If the input set has bounds, the algebraic identity applies whenever the compared inputs are defined; all inputs must also be admissible for a feasible experiment. Bounds do not turn an arbitrary nonlinear actuator into a piecewise linear one. ### Counterfactual Decomposition Where the zero-input branch is defined, $$ f(x,u,t)=f(x,0,t)+G(x,t)u. $$ {#eq-counterfactual} This is a decomposition evaluated at **one state**. Along an actual motion, integrating the two vector-valued terms in a common coordinate chart reconstructs that motion's coordinate increment. The integral of $f(x(t),0,t)$ generally does not equal a forward zero-input motion, which evaluates $f$ at its own evolving state. A coordinate increment, a scalar speed, work and a changed impact outcome are different quantities. ## Time Evolution and Moving Tangent Spaces ### Why Global Superposition Fails Along a reference motion the coefficients are $$ A(t)=f_x(\bar x(t),\bar u(t),t),\qquad B(t)=f_u(\bar x(t),\bar u(t),t). $$ {#eq-time-varying-jacobians} Changing the finite input generally changes the entire state history and hence these coefficients. Integrating a derivative on the old history predicts the first variation of the new history, not its complete finite change. ### Flow Transport Is Not Parallel Transport On $\mathbb R^n$, tangent vectors at different points have a canonical identification as ordinary vectors. Adding them is mathematically meaningful. Nevertheless, for $\dot x=x^2$ on the flat real line, $$ \varphi_T(x_0)=\frac{x_0}{1-Tx_0},\qquad D\varphi_T(x_0)=\frac{1}{(1-Tx_0)^2}, $$ on an interval before the denominator vanishes. The flow is nonlinear although the space has zero intrinsic curvature. Its derivative transports a perturbation according to the dynamic law. Euclidean parallel transport instead leaves vector components constant; it is a different operation. More generally, $T\varphi_{t,s}$ transports initial tangent perturbations to the current state along a fixed prescribed input history. It does not require introducing a metric connection. If a connection is used to compare vectors or define norms, state the choice and do not equate its parallel transport with the dynamic solution derivative. ### Separate the Two Refinement Experiments For sufficiently smooth time dependence, the full nonlinear solution has $$ \begin{aligned} x(t+h)&=x(t)+h f(t)\\ &\quad+\tfrac12h^2\bigl(f_x f+f_u\dot u+f_t\bigr)+O(h^3). \end{aligned} $$ {#eq-taylor-time} Here $f(t)=f(x(t),u(t),t)$. A bounded third time derivative is a sufficient condition for the displayed remainder. $C^1$ state dependence and arbitrary input histories do not provide it. At an input jump, integrate on each time segment with the appropriate one-sided behavior. Euler's method recomputes the **full nonlinear** rate at its numerical state. Under suitable regularity, consistency and stability on a common finite interval, reducing its step size $h$ converges to the nonlinear solution. Integrating the first-variation equation with smaller $h$ instead converges to the first variation. For the exact scalar flow above, the error of predicting a finite initial increment $\epsilon$ by $\varphi_T(x_0)+\epsilon D\varphi_T(x_0)$ is $$ \frac{T\epsilon^2} {(1-Tx_0)^2(1-Tx_0-T\epsilon)}. $$ {#eq-limit-exactness} With $x_0=1$, $T=0.2$ and $\epsilon=0.1$ in nondimensional variables, the error is about $0.00400641$, even with exact integration. It disappears as amplitude tends to zero, not as a solver step tends to zero at fixed amplitude. ### Variational Equation The precise local statement is therefore $$ \dot\eta=A(t)\eta+B(t)v(t),\qquad x_\epsilon(t)=\bar x(t)+\epsilon\eta(t)+o(\epsilon), $$ {#eq-variational} uniformly on a suitable fixed finite interval. A stronger quadratic remainder requires correspondingly stronger regularity and neighborhood bounds. Verify integration convergence and perturbation-amplitude convergence separately. This is the bridge from exact differentiation to a defensible finite prediction. The standard flow and equilibrium results are discussed in [Hunter's ODE notes, Sections 1.3–1.4](https://www.math.ucdavis.edu/~hunter/m207/m207.pdf). The distinction between directional and norm-controlled functional variation is developed in [Liberzon, Section 1.3.2](https://liberzon.csl.illinois.edu/teaching/cvoc.pdf). The counterexamples above make the specific limitations used here explicit. ## Local Optimal Control: LQR and Lyapunov Stability ### Connection to Linear Quadratic Regulation A local derivative model supports controller design when paired with a stated objective and a separate stability argument. #### The LQR Problem Linear Quadratic Regulation (LQR) solves the following problem: **Given:** - Linear system: $\dot{x} = Ax + Bu$ - Quadratic cost: $J = \int_0^\infty (x^T Q x + u^T R u) dt$ - Matrices $Q \succeq 0$ (positive semidefinite), $R \succ 0$ (positive definite) **Find:** Optimal feedback control law minimizing $J$ **Solution:** Linear state feedback $u = -Kx$ where gain $K$ is computed by solving the **algebraic Riccati equation**: $$ A^T P + P A - P B R^{-1} B^T P + Q = 0 $$ {#eq-riccati} with $K = R^{-1} B^T P$. Standard sufficient assumptions for the stabilizing infinite-horizon solution are stabilizability of $(A,B)$ and detectability of $(Q^{1/2},A)$. Select the stabilizing positive-semidefinite solution; the matrix equation alone does not guarantee that an arbitrary root defines an optimal stabilizer. #### What a Nonlinear LQR Design Establishes At an equilibrium, a differentiable nonlinear system has a linear first variation. If the resulting closed-loop Jacobian is Hurwitz, the smooth nonlinear closed loop is locally asymptotically stable. This is a local result with a model-dependent neighborhood. It does not establish that the linear feedback globally minimizes the nonlinear objective. A $C^2$ cost has a second-order Taylor model. Its linear term vanishes at an unconstrained stationary minimum, but the Hessian may be singular or zero: $\ell(x)=x^4$ has no nonzero quadratic leading term. Along a general nominal trajectory, retain cost gradients and state-input cross derivatives rather than assuming each stage is at a minimum. Along a dynamically feasible moving reference, time-varying LQR propagates a value approximation backward from the terminal cost. Solving independent equilibrium Riccati equations at each instant is a different design. The finite-horizon recursion, nonlinear rollout, input constraints, and stability argument must be checked separately. See the [reviewed optimal-control chapter](../The_Geometry_of_Motion/quarto/ch05_optimal_control.html). #### Gain Scheduling Requires Its Own Check Interpolating locally stabilizing gains does not automatically stabilize the scheduled nonlinear system. For a scheduled Lyapunov metric $P(\rho)$, its derivative contains $\sum_i P_{\rho_i}\dot\rho_i$ as well as the closed-loop Jacobian terms. A common Lyapunov function, appropriate scheduling-rate bounds, or a direct trajectory analysis can supply a certificate. Frozen-point eigenvalues alone do not. ### Lyapunov Stability Analysis #### Lyapunov's Indirect Method One of the most powerful tools for analyzing nonlinear stability is **Lyapunov's indirect method**: ::: {.callout-tip} ### Theorem: Lyapunov's Indirect Method For $f$ continuously differentiable on a neighborhood, consider $\dot{x} = f(x)$ with equilibrium $x_e$ (i.e., $f(x_e) = 0$). Let $A = \frac{\partial f}{\partial x}(x_e)$ be the Jacobian at equilibrium. 1. If all eigenvalues of $A$ have strictly negative real parts (Re$(\lambda_i) < 0$), then $x_e$ is **locally asymptotically stable** 2. If any eigenvalue of $A$ has strictly positive real part (Re$(\lambda_j) > 0$), then $x_e$ is **unstable** 3. If all eigenvalues have non-positive real parts with some on imaginary axis, stability is **inconclusive** (higher-order analysis required) ::: #### Why This Works For the stable case, write the displacement dynamics as $\dot z=Az+r(z)$, with $r(z)=o(\|z\|)$. If $A$ is Hurwitz, choose $Q\succ0$ and solve $A^\top P+PA=-Q$ for $P\succ0$. Then $V=z^\top Pz$ satisfies $$ \dot V=-z^\top Qz+2z^\top P r(z). $$ For a sufficiently small neighborhood, the remainder obeys $2\|P\|\|r(z)\|\leq\tfrac12\lambda_{\min}(Q)\|z\|$. Consequently $\dot V\leq-\tfrac12\lambda_{\min}(Q)\|z\|^2$ there. A sufficiently small invariant sublevel set supplies the local stability domain. The proof uses a spectral margin and a bound on the nonlinear remainder; exact differentiation alone is not the certificate. ::: {.callout-note} ### The Pendulum Example Revisited Recall the pendulum at $x_e = [0,0]^T$ has: $$ A = \begin{bmatrix} 0 & 1 \\ -g/\ell & 0 \end{bmatrix}, \quad \lambda = \pm i\sqrt{g/\ell} $$ The eigenvalues are pure imaginary (on the borderline), so Lyapunov's indirect method is **inconclusive**. Indeed, the pendulum equilibrium is: - Stable (nearby trajectories stay nearby) - But not asymptotically stable (they don't converge to equilibrium) - Neutrally stable—consistent with conservative mechanical system For the inverted pendulum at $x_e = [\pi, 0]^T$, one eigenvalue would be positive, confirming instability. ::: #### Design and Validation The pendulum's conserved energy, rather than its imaginary eigenvalues alone, establishes the downward equilibrium's local stability. The inverted equilibrium has a positive eigenvalue and is unstable under the indirect method's hypotheses. For a controller, first check the feasible equilibrium or reference and the controllable directions. Then design for the stated cost, derive a closed-loop neighborhood bound and compare its predicted outputs against nonlinear rollouts. Gain scheduling and actuator saturation require their own checks. A small region in a chosen norm does not guarantee good performance if an output is highly sensitive or the control margin is small. Large excursions, contact transitions and changing shaft states can motivate nonlinear control analysis, but their presence alone does not prove that a particular linear controller will fail. The relevant question is whether the proposed controller meets the declared task tolerances over the required domain. ## Residuals, Reference Motion, and Nonlinearity ### An Interaction Requires Four Comparable Trajectories Fix an initial state, model, time interval, and nominal input $\bar u$. Let $x_{00}$ use $\bar u$, $x_{10}$ use $\bar u+\epsilon v_1$, $x_{01}$ use $\bar u+\epsilon v_2$, and $x_{11}$ use $\bar u+\epsilon(v_1+v_2)$. In a common state coordinate chart, define $$ r(t_1)=x_{11}(t_1)-x_{10}(t_1)-x_{01}(t_1)+x_{00}(t_1). $$ {#eq-thesis-residual-def} The fourth term is the **evolved nominal trajectory**, not the initial state. With zero nominal input it is the evolved zero-input trajectory. For the scalar linear system $\dot x=ax+u$, constant $u$, and duration $T$, the endpoint is $e^{aT}x_0+u(e^{aT}-1)/a$ for $a\ne0$. The four-trajectory residual is zero. Replacing its last term by $x_0$ gives $x_0(1-e^{aT})$, a false interaction even in this linear system. For $\dot x=u^2$ on the Euclidean line and constant perturbations, the same residual is $2\epsilon^2v_1v_2T$. The state space is flat; the input-to-flow map is nonlinear. Ordinary Hessians depend on coordinates and are not intrinsic Riemann curvature tensors. On a manifold, state subtraction requires a declared chart or retraction and a consistent comparison point. ### First Variations and the Mixed Second Variation For a fixed finite horizon, assume $f(x,u,t)$ is $C^2$ in state and input on a neighborhood containing the relevant trajectories, with sufficient time regularity and bounded input directions. Let $$ \begin{aligned} \dot\xi_i&=A(t)\xi_i+B(t)v_i,\qquad \xi_i(t_0)=0,\\ A&=f_x(\bar x,\bar u,t),\qquad B=f_u(\bar x,\bar u,t). \end{aligned} $$ These are exact first variations along one nominal trajectory. Their sum obeys the same linear equation. Define $\eta_{12}$ as the mixed second variation of the solution with respect to the two input amplitudes. It satisfies $$ \dot\eta_{12}=A\eta_{12}+S_{12},\qquad \eta_{12}(t_0)=0, $$ $$ S_{12}=f_{xx}[\xi_1,\xi_2]+f_{xu}[\xi_1,v_2] +f_{xu}[\xi_2,v_1]+f_{uu}[v_1,v_2]. $$ {#eq-curvature-connection} Every derivative is evaluated along the same nominal state and input. For a control-affine system, $f_{uu}=0$, but $f_{xu}$ can be nonzero when the input vector fields depend on state. The interaction is $$ \begin{aligned} r(t_1)&=\epsilon^2\eta_{12}(t_1)+o(\epsilon^2),\\ \eta_{12}(t_1)&=\int_{t_0}^{t_1}\Phi(t_1,t)S_{12}(t)\,dt. \end{aligned} $$ {#eq-residual-scaling} There is no factor of $1/2$ on this mixed interaction: the two symmetric mixed terms in the Taylor expansion cancel that factor. A single trajectory's second-order Taylor error is a different quantity. ### Transport Is Part of the Bound For a chosen consistent norm and fixed interval, $$ \|r(t_1)\|\leq\epsilon^2\int_{t_0}^{t_1} \|\Phi(t_1,t)\|\,\|S_{12}(t)\|\,dt+o(\epsilon^2). $$ {#eq-thesis-residual-bound} The Jacobian matters through $\Phi$ and the first variations inside $S_{12}$. For $\dot x=ax+u^2$, constant directions about $\bar u=0$ give $r=2\epsilon^2v_1v_2(e^{aT}-1)/a$. A constant Hessian source can therefore be amplified exponentially. A bound proportional to $T$ is valid only when its coefficient controls this transport and the other terms over the entire stated interval; it is not a universal linear growth law. ### What a Small Residual Does Not Establish A small interaction residual does not bound every prediction error. For instance, two input directions may have no mutual interaction while each has its own nonlinear self-effect. A zero denominator also makes $\|r\|/\|\delta x\|$ unusable. Choose physical scales or a metric, a reference magnitude bounded away from zero, and a task tolerance: $$ \frac{\|r\|_W}{s_{\mathrm{reference}}}\ll1,\qquad s_{\mathrm{reference}}>0. $$ {#eq-validity-criterion} This is an interaction diagnostic, not a stability or optimality certificate. Separately check finite-model prediction errors, disturbances, discretization error, and active contact modes. A computation should demonstrate convergence under tighter integration tolerances before calling its residual physical. ## Part II: Integral Superposition and Transport {.unnumbered} ## Fundamentals of Integral Accumulation ### Integration Preserves Linearity If $g_1,g_2$ are integrable vector-valued functions in one coordinate representation and $a,b$ are fixed scalars, then $\int(ag_1+bg_2)=a\int g_1+b\int g_2$. No dynamic linearity assumption is needed for this algebraic property. The substantive question is whether the integrands being added describe the same experiment. For first variations along one common reference, $$ \dot\eta=A(t)\eta+B(t)v(t),\qquad \eta(t_0)=\eta_0. $$ {#eq-variational-repeat} The coefficients are fixed by the reference. The derivative direction $v$ may vary over time. Superposition holds between solutions of this linear variational system when initial variations and input directions are combined consistently. ### The Integration Theorem Direct integration gives the exact identity $$ \eta(T)-\eta(t_0) =\int_{t_0}^{T}\bigl(A(t)\eta(t)+B(t)v(t)\bigr)\,dt. $$ {#eq-integrated-variation} The first term includes the unknown evolving variation; omitting it is generally incorrect. Introduce the state-transition matrix $\Phi(t,s)$ by $\partial_t\Phi(t,s)=A(t)\Phi(t,s)$ and $\Phi(s,s)=I$. Variation of constants gives $$ \eta(T)=\Phi(T,t_0)\eta_0+ \int_{t_0}^{T}\Phi(T,t)B(t)v(t)\,dt. $$ {#eq-endpoint-variation} To verify the formula, differentiate its right-hand side with respect to $T$. The integral's moving upper limit supplies $B(T)v(T)$, and the derivative of each transition matrix supplies $A(T)$ times the complete expression. Its initial value is $\eta_0$. Uniqueness of the linear initial-value problem then identifies it with $\eta$. ### What Superposes: Variations About One Reference Let $x_{\epsilon,i}$ use input $\bar u+\epsilon v_i$ from the common initial state, and let $x_{00}$ use $\bar u$. The first variation is $$ \eta_i(T)=\left.\partial_\epsilon x_{\epsilon,i}(T)\right|_0, \qquad x_{\epsilon,i}(T)-x_{00}(T)=\epsilon\eta_i(T)+o(\epsilon). $$ It is not $x_{\epsilon,i}(T)-x(t_0)$. Subtracting the initial state mixes the nominal displacement with the perturbation response. The sum $\eta_1+\eta_2$ is exactly the derivative for $v_1+v_2$; the finite branched responses retain their nonlinear residual. When an input or initial condition is uncertain, propagate its joint uncertainty through the common derivative. Summing separate output variances is justified only when the required cross-covariances vanish. Linearity of the variational equation does not establish independence of measured joint torques. ## State Transition Operator ### Solution of the Variational Equation For zero input variation, $$ \eta(t)=\Phi(t,s)\eta(s). $$ {#eq-state-transition-solution} For finite-dimensional ODEs with continuous coefficients, $\Phi$ is invertible and satisfies $\Phi(t,r)\Phi(r,s)=\Phi(t,s)$ and $\Phi(s,t)=\Phi(t,s)^{-1}$. This is a two-time composition, or cocycle, property. A time-varying system generally depends on both endpoint times, not just their difference. ### Interpretation: Transport Then Accumulate An input direction applied near time $t$ changes the rate through $B(t)$. The resulting state variation then evolves through $\Phi(T,t)$. A small instantaneous gain can therefore be followed by strong amplification; a large acceleration increment can be followed by cancellation in the selected endpoint output. Transition matrices at different times need not commute. Take nondimensional matrices $$ A_1=\begin{pmatrix}0&1\\0&0\end{pmatrix}, \qquad A_2=\begin{pmatrix}0&0\\1&0\end{pmatrix}, $$ on two consecutive unit-duration intervals. Their ordered transition is $$ e^{A_2}e^{A_1}=\begin{pmatrix}1&1\\1&2\end{pmatrix}, $$ whereas reversing their order gives $\begin{pmatrix}2&1\\1&1\end{pmatrix}$. Neither equals $e^{A_1+A_2}$. Replacing $\Phi$ by an exponential of an integrated Jacobian needs additional commutation conditions. ### Time-Invariant Special Case For constant $A$, $$ \Phi(t,s)=e^{A(t-s)}. $$ {#eq-matrix-exponential} The LTI solution is $$ x(T)=e^{A(T-t_0)}x_0+\int_{t_0}^{T}e^{A(T-t)}Bu(t)\,dt. $$ {#eq-lti-solution} This is the complete motion for the linear plant. The same formula with time-varying reference coefficients gives the derivative response for a nonlinear plant; the identical notation does not make the latter a complete nonlinear solution. ## Variational Dynamics and Sensitivity Analysis ### Sensitivity to Initial Conditions Holding the prescribed input history fixed, define $$ S_x(t)=D_{x_0}x(t;x_0,\bar u). $$ {#eq-sensitivity-ic} Then $$ S_x(t)=\Phi(t,t_0), $$ {#eq-sensitivity-phi} $$ \dot S_x=A(t)S_x,\qquad S_x(t_0)=I. $$ {#eq-sensitivity-evolution} For a differentiable output $y=h(x(T))$, its initial-state derivative is $Dh(\bar x(T))\Phi(T,t_0)$. This task map is essential when the state contains angles, velocities and shaft coordinates while the question concerns clubface orientation or impact position. ### An Input Kernel Is a Density At fixed terminal time and initial state, write $$ K(T,t)=\Phi(T,t)B(t). $$ {#eq-sensitivity-input} The kernel is a **sensitivity density**, to be integrated against an input perturbation and time. If input has torque units, its entries have corresponding state-per-torque-per-time units. It is not generally the derivative of the endpoint with respect to changing a finite torque value at one isolated instant. For bounded measurable input histories, changing the input at a single time leaves its equivalence class almost everywhere unchanged and does not change the integral solution. A short rectangular pulse of height $\epsilon a$ on $[\tau,\tau+h]$ instead produces the first-order endpoint change $\epsilon\int_\tau^{\tau+h}K(T,t)a\,dt$. For continuous $K$ this approaches $\epsilon h K(T,\tau)a$ as $h$ tends to zero. For the nondimensional system $\dot x=-x+u$, $x(0)=0$, the exact pulse response is $$ \Delta x(T)=\epsilon a\,e^{-(T-\tau)}(e^h-1). $$ At $T=1$, $\tau=0.3$, $\epsilon a=0.2$ and $h=10^{-4}$, it is about $9.9322\times10^{-6}$, not the height-only value $0.2e^{-0.7}$. A Dirac impulse has a specified area and different units; finite impulsive interventions can leave the small-perturbation regime and need their own dynamic interpretation. ### The Endpoint Derivative Is an Operator Let $\mathcal E(u)=x(T;x_0,u)$ for a fixed interval and initial state. Under uniform smoothness and neighborhood bounds that make this solution map Fréchet differentiable for bounded input histories, its derivative acts on a whole function: $$ D\mathcal E(\bar u)[v] =\int_{t_0}^{T}\Phi(T,t)B(t)v(t)\,dt. $$ {#eq-integrated-sensitivity} Writing only $\int\Phi B\,dt$ omits the argument. That matrix is the response to a **constant** input-offset parameter, not the entire functional derivative. In the scalar example at $T=1$, $v(t)=1$ gives $1-e^{-1}\approx0.6321$, while $v(t)=t$ gives $e^{-1}\approx0.3679$. Timing changes the endpoint effect even though both directions use the same kernel. This is a precise reason that the sequence of torque and posture changes matters in the swing. ## Adjoint (Costate) Dynamics ### Backward Sensitivity For fixed $t_0,T,x_0$, consider $$ J(u)=\phi(x(T))+\int_{t_0}^{T}L(x(t),u(t),t)\,dt. $$ {#eq-cost-functional} Here $L$ is the chosen running cost, not necessarily a mechanical Lagrangian. Use column gradients. Along the prescribed nominal history define an adjoint $\lambda$ so that $$ \begin{aligned} C(t,x)&=\phi(X(T))+\int_t^T L(X(s),\bar u(s),s)\,ds,\\ \lambda(t)&=\left.\nabla_x C(t,x)\right|_{x=\bar x(t)}. \end{aligned} $$ {#eq-costate-def} Here $X(s)=x(s;t,x,\bar u)$ starts from state $x$ at time $t$. This is the derivative of the future cost when the **future input history is held fixed**. It is not automatically the gradient of the optimally minimized value function. ### Adjoint Evolution Equation The adjoint solves $$ \dot\lambda=-A^\top\lambda-L_x,\qquad $$ {#eq-adjoint-evolution} $$ \lambda(T)=\phi_x(\bar x(T)). $$ {#eq-adjoint-terminal} The first cost variation is $$ \delta J=\phi_x^\top\eta(T)+ \int_{t_0}^T(L_x^\top\eta+L_u^\top v)\,dt. $$ Using the forward and backward equations, $\frac{d}{dt}(\lambda^\top\eta)=-L_x^\top\eta+\lambda^\top Bv$. Integrating this identity cancels the state-variation terms and yields $$ DJ(\bar u)[v]=\int_{t_0}^{T} \bigl(L_u+B^\top\lambda\bigr)^\top v\,dt. $$ {#eq-optimality-condition} The quantity $L_u+B^\top\lambda$ is a gradient density for the specified pairing. If the initial state varies, add $\lambda(t_0)^\top\eta_0$. Different inner products on the space of input functions represent the same derivative by different gradients. For an interior unconstrained stationary input, this density vanishes almost everywhere. With bounds or other constraints, admissible directions are restricted and the appropriate variational inequality or optimal-control conditions replace unrestricted stationarity. Endpoint constraints and variable terminal time also alter boundary conditions. ### The Policy Determines the Sensitivity If a controller recomputes $u=\kappa(x,t)$ after a state perturbation, the closed-loop state derivative is $A+B\kappa_x$. The running cost becomes $L(x,\kappa(x,t),t)$ and its state gradient includes $\kappa_x^\top L_u$. For $\dot x=x+u$ with $u=-2x$, the closed-loop perturbation decays as $e^{-t}$; freezing the nominal input history instead gives the open-loop sensitivity $e^t$. Both are valid answers to different experiments. For a golfer, a prescribed torque replay, a model of feedback and a reoptimized motion likewise define different comparisons. An adjoint cannot decide which one represents the player's response. ### A Reproducible Gradient Check Use the nondimensional system $\dot x=x^2+u$, $x(0)=0.2$, $T=0.3$, nominal $u=0.4$, $\phi=x^2/2$ and $L=(x^2+0.1u^2)/2$. Then $\dot\lambda=-2x\lambda-x$, $\lambda(T)=x(T)$, and the gradient density is $\lambda+0.1u$. For direction $v(t)=1+t$, compare its integral against $[J(\bar u+\epsilon v)-J(\bar u-\epsilon v)]/(2\epsilon)$ from two independent nonlinear integrations. The supplied regression checks agreement as well as solver accuracy. This verifies the derivative computation for this example; it does not validate a golf model. ### Task and Event Interpretation An endpoint state sensitivity becomes a clubhead sensitivity only after the appropriate kinematic output derivative is applied. At first contact, the event time usually varies too. For a transversal event $c(x,t)=0$, the time correction has denominator $c_x f+c_t$; a small or zero denominator signals grazing and possible loss of differentiability. Include the event-time derivative, and any impact/reset derivative, before interpreting a pre-impact fixed-time kernel as a ball-strike prediction. The adjoint weights mechanical effects by a declared objective, while the forward model transports them through changing posture and shaft state. Neither produces a unique allocation of observed performance to muscles, segments or coaching intentions. That inference requires a separately defined intervention and evidence. ## Physical Interpretations: Impulse and Work ## Force-Momentum: Impulse Integral ### Newton's Second Law For a mechanical Lagrangian $L_m=T-V$ with $p=\partial L_m/\partial\dot q=M(q)\dot q$, the Euler–Lagrange momentum equation can be written: $$ \frac{dp}{dt} = f(q, \dot{q}) + u $$ {#eq-newton} Here $u$ is the chosen generalized applied load and $f=\partial L_m/\partial q+Q_{\mathrm{other}}$. The derivative $\partial L_m/\partial q$ includes configuration-dependent kinetic terms as well as potential forces. Canonical generalized momentum need not equal a Cartesian momentum component. For Cartesian translation with constant mass, the equation reduces to Newton's law with all other external forces included in $f$. ### Impulse Definition The **impulse** delivered by a force $u$ over $[t_0, t_1]$ is: $$ J = \int_{t_0}^{t_1} u(\tau) \, d\tau $$ {#eq-impulse} From Newton's second law: $$ \int_{t_0}^{t_1} \frac{dp}{dt} dt = p(t_1) - p(t_0) = \int_{t_0}^{t_1} [f + u] dt $$ {#eq-impulse-momentum} If the system starts at rest ($\dot{q}(t_0) = 0$) and the integral of every retained term in $f$ is demonstrably negligible on the stated interval, then approximately: $$ p(t_1) \approx \int_{t_0}^{t_1} u(\tau) d\tau = J $$ {#eq-impulse-approx} The integral of the complete right-hand side equals the canonical momentum change. The impulse of the selected input alone equals that change only when all remaining integrated terms vanish; in Cartesian coordinates the total external impulse equals the change in linear momentum. ### Superposition of Impulses If two forces $u_1$ and $u_2$ are applied: $$ \begin{aligned} J_{12}&=\int_{t_0}^{t_1}[u_1(\tau)+u_2(\tau)]\,d\tau\\ &=\int_{t_0}^{t_1}u_1(\tau)\,d\tau+ \int_{t_0}^{t_1}u_2(\tau)\,d\tau=J_1+J_2. \end{aligned} $$ {#eq-impulse-superposition} This is **exact**—impulses superpose because integration is linear. The resulting momentum change is: $$ p(t_1) = p(t_0) + J_1 + J_2 + \int_{t_0}^{t_1}f(q(t),\dot q(t))\,dt $$ {#eq-momentum-superposition} This is exact bookkeeping on the combined trajectory. The remaining integral is not generally higher order or small; gravity can contribute at first order in elapsed time. If force laws or feedback depend on the state, their actual histories also change across separate counterfactual runs. ::: {.callout-note} ## Interpretation The complete momentum balance adds retained load terms in one coordinate representation along the actual motion. Its individual input integrals are signed accounting terms, not independent forward-trajectory predictions. ::: ## Power-Energy: Work Integral ### Power and Work For a mechanical system, **power** delivered by a force is: $$ P(t) = u(t)^T \dot{q}(t) $$ {#eq-power} The **work** done over $[t_0, t_1]$ is: $$ W = \int_{t_0}^{t_1} P(\tau) d\tau = \int_{t_0}^{t_1} u(\tau)^T \dot{q}(\tau) d\tau $$ {#eq-work} For generalized load $u$ paired with $\dot q$, this is the work of that load. The kinetic-energy balance must include the work $W_{\mathrm{other}}$ of all other forces, including gravity when it is not absorbed into potential energy: $$ W+W_{\mathrm{other}}=E_k(t_1)-E_k(t_0) $$ {#eq-work-energy} Here $E_k=\tfrac12\dot q^TM\dot q$. For a time-independent conservative potential and only dissipative loads $-d$ besides $u$, the mechanical-energy balance is $\Delta(E_k+V)=\int(u-d)^T\dot q\,dt$. Ideal stationary constraint reactions do no work; moving supports can. These assumptions must be stated before equating actuator work with an energy change. ### Superposition of Work If two forces $u_1$ and $u_2$ are applied: $$ \begin{aligned} W_{12}&=\int_{t_0}^{t_1}[u_1+u_2]^\top\dot q\,d\tau\\ &=\int_{t_0}^{t_1}u_1^\top\dot q\,d\tau+ \int_{t_0}^{t_1}u_2^\top\dot q\,d\tau=W_1+W_2. \end{aligned} $$ {#eq-work-superposition} Again, **exact** additivity due to integral linearity. The velocity history is the same combined-motion history in both work terms. For a unit mass initially at rest, two constant unit forces over a unit interval produce velocity $2t$, final kinetic energy 2 and work 1 from each force. Running either force alone gives velocity $t$ and work $1/2$. The separate-run works sum to 1, not 2. This interaction occurs even with linear force-to-velocity dynamics because power multiplies force by velocity. ::: {.callout-warning} ## Variation vs. Causation The force-work terms add along the same combined trajectory, providing a valid signed energy budget when all terms are included. Those terms are not the works obtained by running each input alone, and normalized work fractions do not uniquely identify physiological causation or efficiency. Negative work and an initial energy store can also make naive percentage interpretations misleading. Work accounting and first-variation transport are distinct linear operations. Neither automatically assigns a unique fraction of an observed outcome to a physical cause. ::: ## General Principle: Integral Accumulation Both impulse and work illustrate a general principle: ::: {.callout-important} ## Integral Accumulation Principle **Along a fixed trajectory** $x(t)$, for any quantity defined as an integral of a linear function of inputs: $$ \Psi = \int_{t_0}^{t_1} g(x(t), u(t), t) \, dt $$ where $g$ is linear in $u$ with $x(t)$ held fixed, the contributions from different inputs superpose exactly: $$ \Psi(u_1 + u_2) = \Psi(u_1) + \Psi(u_2) $$ **Important caveat:** This integral identity holds for any such input histories evaluated along the same fixed trajectory $x(t)$; the histories need not be infinitesimal. If $x(t)$ itself depends nonlinearly on $u$ (as it generally does in closed-loop or trajectory optimization settings), then the principle does not hold globally — changing $u$ changes the trajectory, and the two contributions cannot be evaluated independently. ::: Impulse and work are special cases where $g = u$ and $g = u^T \dot{q}$ respectively. The linearity of integration ensures exact additivity of these quantities, even as the underlying dynamics remain nonlinear. \newpage ## Finite-Horizon Residuals and Error Quantification ### Keep the Baseline and the Perturbation Definition Fixed The four-trajectory interaction in @eq-thesis-residual-def can also be written using finite deviations $\Delta x_{ij}=x_{ij}-x_{00}$: $$ r(t_1)=x_{11}(t_1)-x_{10}(t_1)-x_{01}(t_1)+x_{00}(t_1). $$ {#eq-residual-global} $$ r=\Delta x_{11}-\Delta x_{10}-\Delta x_{01}. $$ {#eq-residual-variation} These finite deviations differ from the first variations $\xi_i$, whose equation is linear and whose interaction is exactly zero. The notation must not quietly turn one into the other. For a smooth solution map on a fixed interval, $$ \|r\|=O(\epsilon^2). $$ {#eq-residual-order} If $\|\Phi(t_1,t)\|\,\|S_{12}(t)\|\leq C_T$ throughout an interval of duration $T$, then $$ \|r\|\leq\epsilon^2 C_T T+o(\epsilon^2). $$ {#eq-residual-bound-time} $C_T$ can grow rapidly with the horizon. Quadratic order does not mean exact quadrupling at arbitrary amplitudes: the mixed coefficient can vanish, higher orders can matter, and contacts can destroy smooth dependence. Convergence in perturbation amplitude should be tested at a fixed horizon. ### Interaction and Single-Trajectory Approximation Error The relevant Hessians are $$ H_x=f_{xx},\qquad H_u=f_{uu},\qquad H_{xu}=f_{xu}. $$ {#eq-hessians} The mixed interaction uses $S_{12}$ from @eq-curvature-connection: $$ r(t_1)=\epsilon^2\int_{t_0}^{t_1}\Phi(t_1,t)S_{12}(t)\,dt+o(\epsilon^2). $$ {#eq-residual-hessian} By comparison, a single-input expansion $x_\epsilon=\bar x+\epsilon\xi+\tfrac12\epsilon^2\eta+o(\epsilon^2)$ has a self-variation source $f_{xx}[\xi,\xi]+2f_{xu}[\xi,v]+f_{uu}[v,v]$. Both calculations need transport. Neither identifies intrinsic manifold curvature from a measured residual. ### A Task-Based Interpretation for the Swing Specify whether comparisons use the same elapsed time or an event such as first ball contact. A perturbed trajectory may reach that event at a different time, so a contact-event endpoint map must include event-time sensitivity; smooth-flow expressions alone do not provide it. A changed contact mode requires a hybrid analysis. In a common coordinate chart, compare residuals using declared scales: $$ \chi=\|r\|_W/s_{\mathrm{reference}},\qquad s_{\mathrm{reference}}>0. $$ {#eq-validity-quantitative} If the chosen perturbation response scale is proportional to $\epsilon$, a bounded mixed response gives a diagnostic of order $$ \chi=O(\epsilon C_T T). $$ {#eq-validity-simplified} For a prescribed tolerance $\chi_{\max}$, a leading-order estimate is $$ \epsilon C_T T\lesssim\chi_{\max}, $$ {#eq-validity-epsilon} with scales absorbed into $C_T$. It is meaningful only where the derivative bounds and remainder estimates hold. There is no universal threshold such as 0.1 at which a nonlinear controller becomes necessary. Applied to a swing model, this procedure can quantify how two chosen input interventions interact in club position, face orientation, or speed. It does not recover a unique causal allocation of observed performance to muscles. Geometry, state-dependent leverage, activation, shaft storage, feedback, and impact timing can all alter the endpoint map. Compare those mechanisms through explicit model variants and measured uncertainties. ## Part III: Hamiltonian Structure and Optimal Control {.unnumbered} ## Variational Principles and Hamiltonian Formulation ### Mechanics and Optimization Ask Different Questions Hamilton's stationary-action principle determines conservative mechanical trajectories by varying configurations with fixed endpoints in the action $\int L_m(q,\dot q,t)\,dt$. The stationary action need not be a minimum. Applied nonconservative forces enter through the Lagrange–d'Alembert formulation. Optimal control instead chooses an admissible input to minimize a specified objective subject to dynamics. Its running cost need not be kinetic minus potential energy, and the mechanical system does not generally minimize it without an appropriately designed input. Here use $\ell$ for that cost and $L_m$ for a mechanical Lagrangian. For a fixed initial state, the objective is $$ J(u)=\phi(x(T),T)+\int_{t_0}^T\ell(x,u,t)\,dt. $$ {#eq-action} The state equation is $\dot x=f(x,u,t)$, and admissible inputs, endpoint conditions and path constraints must be specified. All gradients below are column vectors. ### A Consistent Minimization Convention For a normal extremal, define $$ H(x,u,\lambda,t)=\ell(x,u,t)+\lambda^\top f(x,u,t). $$ {#eq-hamiltonian} The corresponding necessary conditions are $$ \dot x=H_\lambda=f,\qquad \dot\lambda=-H_x,\qquad u^*(t)\in\underset{u\in\mathcal U}{\operatorname{argmin}}H. $$ {#eq-hamiltonian-eqs} The pointwise condition holds almost everywhere under the hypotheses of the relevant optimal-control principle. For an unconstrained interior solution, it implies $H_u=0$. Stationarity is weaker than minimization, and necessary conditions alone do not prove optimality. Some texts use a Hamiltonian with **negative** running cost and maximize it. For example, [Liberzon, Section 3.4](https://liberzon.csl.illinois.edu/teaching/cvoc.pdf), uses $p^\top f-\ell$. Setting $p=-\lambda$ makes that Hamiltonian the negative of ours. Both conventions are consistent; combining our positive-cost Hamiltonian with that maximization rule reverses the result. The general multiplier formulation includes a nonnegative cost multiplier $\lambda_0$ in $H=\lambda_0\ell+\lambda^\top f$, with nontrivial multipliers. The normalization $\lambda_0=1$ applies to normal extremals. Abnormal extremals with $\lambda_0=0$ require separate treatment; they are not excluded merely by writing the normal equations. ### Costate as Sensitivity Along a fixed nominal input history, the adjoint represents the gradient of the future cost defined in Part II. At an optimum, it can also agree with the gradient of a differentiable value function when the associated dynamic-programming conditions hold. Neither interpretation makes a general costate a physical momentum. For a fixed terminal time and a free terminal state with only terminal cost, $$ \lambda(T)=\phi_x(x(T),T). $$ {#eq-terminal-costate} Backward integration also needs the nominal state and input histories: a terminal state and an objective alone do not supply all coefficients of the adjoint equation. ### Deriving the Dynamic Residual Correctly Consider the fixed-time, unconstrained-interior problem $$ \min_u\ \phi(x(T),T)+\int_{t_0}^T\ell(x,u,t)\,dt $$ {#eq-optimal-control-problem} subject to $$ \dot x=f(x,u,t),\qquad x(t_0)=x_0. $$ {#eq-dynamics-constraint} Introduce a multiplier for the complete dynamic residual: $$ \widetilde J=\phi(x(T),T)+ \int_{t_0}^T\left[\ell+\lambda^\top(f-\dot x)\right]\,dt. $$ {#eq-augmented-functional} With $$ H=\ell+\lambda^\top f, $$ {#eq-hamiltonian-def} this is $$ \widetilde J=\phi+\int_{t_0}^T(H-\lambda^\top\dot x)\,dt. $$ {#eq-augmented-hamiltonian} Vary $x,u,\lambda$ independently, with $\delta x(t_0)=0$. The first variation is $$ \begin{aligned} \delta\widetilde J={}&\phi_x^\top\delta x(T)\\ &+\int_{t_0}^T\left[ H_x^\top\delta x+H_u^\top\delta u +(f-\dot x)^\top\delta\lambda -\lambda^\top\delta\dot x\right]dt. \end{aligned} $$ {#eq-first-variation} The coefficient of $\delta\lambda$ is **$f-\dot x$**, not just $H_\lambda=f$. Omitting the second term would impose $f=0$ instead of the state equation. Integration by parts gives $$ -\int_{t_0}^T\lambda^\top\delta\dot x\,dt =-\left[\lambda^\top\delta x\right]_{t_0}^{T} +\int_{t_0}^T\dot\lambda^\top\delta x\,dt. $$ {#eq-integration-by-parts} Consequently, $$ \begin{aligned} \delta\widetilde J={}&(\phi_x-\lambda(T))^\top\delta x(T)\\ &+\int_{t_0}^T\left[ (H_x+\dot\lambda)^\top\delta x +H_u^\top\delta u +(f-\dot x)^\top\delta\lambda\right]dt. \end{aligned} $$ {#eq-combined-variation} Vanishing coefficients yield the state equation, costate equation, interior stationarity and free-state terminal condition. A fixed terminal state has $\delta x(T)=0$, so its multiplier is not determined by that free-state condition. Bounded controls restrict the admissible variations; needle variations underpin the stronger pointwise minimum principle. ### Endpoint and Time Conditions With a regular endpoint equality $\psi(x(T),T)=0$, an endpoint multiplier $\nu$ gives $$ \lambda(T)=\phi_x+\psi_x^\top\nu. $$ When terminal time is also free and the usual regularity conditions apply, $$ H(T)+\phi_T+\nu^\top\psi_T=0. $$ These are partial derivatives of the explicit endpoint functions. For an autonomous minimum-time problem with fixed terminal state and no terminal cost, the normal Hamiltonian is zero along an extremal. A fixed terminal time does not supply this zero-Hamiltonian condition. ### Bounded Inputs and Bang-Bang Control For our minimization convention, $$ u^*(t)\in\underset{u\in\mathcal U}{\operatorname{argmin}} H(x^*,u,\lambda^*,t). $$ {#eq-pontryagin-maximum} If $H$ is affine in a scalar $u\in[-1,1]$, a nonzero coefficient selects an endpoint. If the coefficient vanishes over an interval, the minimum condition alone does not select the input; a singular arc may require higher-order conditions. Thus bounded inputs do not universally imply a specified bang-bang sequence. ### Minimum Squared Input Is Not Mechanical Energy For the LTI system $\dot x=Ax+Bu$, fixed endpoints and $R\succ0$, consider $$ J=\tfrac12\int_{t_0}^T u^\top Ru\,dt. $$ {#eq-min-energy-cost} This is a weighted squared-input objective. Calling it energy requires a separate physical model relating the input to power. The Hamiltonian and stationary input are $$ H=\tfrac12u^\top Ru+\lambda^\top(Ax+Bu), $$ {#eq-min-energy-hamiltonian} $$ u^*=-R^{-1}B^\top\lambda, $$ {#eq-min-energy-control} $$ \dot\lambda=-A^\top\lambda. $$ {#eq-min-energy-costate} The state/costate boundary-value system is $$ \frac{d}{dt}\begin{pmatrix}x\\\lambda\end{pmatrix} =\begin{pmatrix}A&-BR^{-1}B^\top\\0&-A^\top\end{pmatrix} \begin{pmatrix}x\\\lambda\end{pmatrix}. $$ {#eq-min-energy-tpbvp} Let $d=x_T-e^{A(T-t_0)}x_0$ and $W=\int_{t_0}^T e^{A(T-t)}BR^{-1}B^\top e^{A^\top(T-t)}\,dt$. If $W$ is nonsingular, the optimal input is $u^*(t)=R^{-1}B^\top e^{A^\top(T-t)}W^{-1}d$ and $J^*=\tfrac12d^\top W^{-1}d$. For singular $W$, feasibility requires $d$ in its range; the corresponding minimum uses the pseudoinverse on that reachable subspace. For the nondimensional double integrator from rest to rest over displacement $D$ in fixed time $T$, this gives $u(t)=6D/T^2-12Dt/T^3$ and $J=6D^2/T^3$ when $R=1$. The checked example $D=3,T=2$ has zero final velocity, displacement 3 and cost 6.75. If input bounds exclude that history, it is not the constrained optimum. ### A Complete Minimum-Time Example For $\dot q=v$, $\dot v=u$, $|u|\leq1$, transfer $(q,v)=(0,0)$ to $(1,0)$ in minimum nondimensional time: $$ J=T=\int_0^T1\,dt, $$ {#eq-time-optimal-cost} $$ H=1+\lambda_qv+\lambda_vu. $$ {#eq-time-optimal-hamiltonian} The minimizing rule is $$ u^*=-\operatorname{sign}(\lambda_v),\qquad\lambda_v\ne0. $$ {#eq-time-optimal-control} The solution uses $u=1$ for $0<t<1$ and $u=-1$ for $1<t<2$. Its velocity is $t$ then $2-t$, so the traveled displacement is 1. For any rest-to-rest history of duration $T$, $v(t)\leq\min(t,T-t)$ and hence $q(T)\leq T^2/4$. The achieved duration 2 therefore meets a global lower bound for this specific problem. Choose $\lambda_q=-1$, $\lambda_v=t-1$. Then $$ H=1-v+(t-1)u=0 $$ {#eq-double-integrator-hamiltonian} on both arcs, verifying the normal free-time condition. The input value at the isolated switching instant does not change the trajectory. This example's single switch is not a theorem about arbitrary nonlinear plants. ### LQR as a Special Boundary-Value Problem For $$ J=\tfrac12x(T)^\top Q_f x(T) +\tfrac12\int_{t_0}^T(x^\top Qx+u^\top Ru)\,dt, $$ {#eq-lqr-cost} $$ \dot x=Ax+Bu, $$ {#eq-lqr-dynamics} the Hamiltonian is $$ H=\tfrac12(x^\top Qx+u^\top Ru)+\lambda^\top(Ax+Bu). $$ {#eq-lqr-hamiltonian} With $Q,Q_f\succeq0$ and $R\succ0$, $$ u^*=-R^{-1}B^\top\lambda, $$ {#eq-lqr-optimal-control} $$ \dot\lambda=-Qx-A^\top\lambda. $$ {#eq-lqr-costate-eq} Substituting $\lambda=P(t)x$ and matching coefficients yields $$ -\dot P=Q+A^\top P+PA-PBR^{-1}B^\top P, \qquad P(T)=Q_f, $$ {#eq-riccati-de} $$ u^*=-R^{-1}B^\top P(t)x. $$ {#eq-lqr-feedback} This solves the stated finite-horizon linear-quadratic problem. For a nonlinear plant it supplies a local approximation or controller design, whose finite performance and stability require separate verification. ### Box Constraints Require the Correct Metric For $u_{\min}\leq u\leq u_{\max}$, use inequalities $u_{\min}-u\leq0$ and $u-u_{\max}\leq0$. The stationarity condition is $$ H_u-\mu_{\min}+\mu_{\max}=0, $$ {#eq-kkt-optimality} with nonnegative multipliers and componentwise complementarity, $$ \mu_{\min}\odot(u-u_{\min})=0,\qquad \mu_{\max}\odot(u_{\max}-u)=0. $$ {#eq-thesis-complementarity} For a convex quadratic Hamiltonian, minimization is projection of the unconstrained solution in the quadratic weight's metric. Componentwise clipping solves this projection for a positive diagonal weight and box constraints; it generally fails for a coupled weight. For example, minimize $\tfrac12(u-a)^\top R(u-a)$ with $a=(2,0)$, $R=\begin{pmatrix}2&1\\1&2\end{pmatrix}$ and $u\in[-1,1]^2$. Clipping gives $(1,0)$, while the actual optimum is $(1,1/2)$. Its gradient is $(-3/2,0)$, balanced by upper-bound multiplier $(3/2,0)$. Changing one saturated component changes the best value of the other through the cross term. ### Path Constraints and Their Limits For regular mixed state/input inequalities $c(x,u,t)\leq0$, a smooth active-arc multiplier formulation can use $$ H_{\mathrm{aug}}=H+\mu^\top c, $$ {#eq-augmented-ham-state} $$ \dot\lambda=-H_x-c_x^\top\mu, $$ {#eq-costate-with-constraints} together with stationarity, feasibility, nonnegative multipliers and complementarity. Applicable constraint qualifications are required. Pure state constraints can involve boundary arcs, differentiated tangency conditions and costate jumps or measure multipliers. They are not universally handled by appending an ordinary smooth multiplier to these equations. A numerical transcription must check constraints between nodes as well as at nodes, and a contact reset needs its own model. These distinctions matter whenever a swing model imposes grip closure, ground contact, actuator bounds or a precise impact event. ## Numerical Implementation and Discrete-Time Systems ### Differentiate the Actual Discrete Transition Let $F_k(x,u)$ be the transition used by the implementation, with an explicitly specified input hold over $[t_k,t_{k+1}]$. It may be the exact held-input flow or a numerical approximation. A nominal rollout satisfies $\bar x_{k+1}=F_k(\bar x_k,\bar u_k)$. Euler defines $$ F_k^E(x,u)=x+h_kf(x,u,t_k). $$ {#eq-euler} The expansion about that nominal transition is $$ \begin{aligned} x_{k+1}={}&\bar x_{k+1}+A_k^d\delta x_k+B_k^d\delta u_k\\ &+o(\|(\delta x_k,\delta u_k)\|). \end{aligned} $$ {#eq-discrete-linearization} Here $$ A_k^d=D_xF_k,\qquad B_k^d=D_uF_k. $$ {#eq-discrete-jacobians} For Euler these are $I+h_kf_x$ and $h_kf_u$, respectively. For an RK4 rollout, differentiate its entire transition, including all intermediate evaluations. Using Euler derivatives with an RK4 rollout is an additional approximation, not the exact derivative of the implemented model. An arbitrarily drawn state path has defects $F_k(\bar x_k,\bar u_k)-\bar x_{k+1}$; either include those defects in the optimization or generate the nominal states by rollout. The first variation of the chosen discrete model obeys exactly $$ \eta_{k+1}=A_k^d\eta_k+B_k^dv_k. $$ {#eq-discrete-variation} This exact derivative statement does not assert exact finite perturbation prediction. Even a continuous system affine in held input can have a nonlinear finite-step map in that input. ### Discrete Transport and Input Holding Define $$ \Phi_{k+1|k}=A_k^d, $$ {#eq-discrete-phi} $$ \Phi_{N|0}=A_{N-1}^d\cdots A_1^dA_0^d. $$ {#eq-phi-chain} Then $\eta_N=\Phi_{N|0}\eta_0+\sum_{k=0}^{N-1}\Phi_{N|k+1}B_k^dv_k$. The empty product is the identity. The order matters. For $\dot q=v,\ \dot v=u$ with constant input over duration $h$, the exact matrices are $A^d=\begin{pmatrix}1&h\\0&1\end{pmatrix}$ and $B^d=(h^2/2,h)^\top$. Euler instead gives $(0,h)^\top$ for the input column. The difference is real displacement accumulated while velocity changes; it is not a change in the dimension of a tangent space. ### Integration Order and Regularity For autonomous dynamics with a held input, classical RK4 uses $$ \begin{aligned} k_1&=f(x_k,u_k),\\ k_2&=f(x_k+\tfrac h2k_1,u_k),\\ k_3&=f(x_k+\tfrac h2k_2,u_k),\\ k_4&=f(x_k+hk_3,u_k),\\ F_k^{RK4}&=x_k+\tfrac h6(k_1+2k_2+2k_3+k_4). \end{aligned} $$ {#eq-rk4} For explicitly time-dependent dynamics, use the corresponding stage times. Input interpolation must also match the specified hold. On a fixed smooth interval with sufficient bounded derivatives and stable error propagation, Euler has one-step state error $O(h^2)$ and global error $O(h)$; RK4 has one-step error $O(h^5)$ and global error $O(h^4)$. Some references divide the one-step defect by $h$ when defining local truncation error; the convention must be stated. An input jump or impact can violate the smoothness needed by a high-order formula. Align integration segments with known discontinuities and locate events to a stated tolerance. Increasing formal order does not repair an incorrect contact model. ### Adaptive Integration Is Not a Stiffness Cure An embedded pair estimates local error from two formulas. A simplified controller for an estimate proportional to $h^5$ is $$ h_{\rm new}=s\,h_{\rm old} \left(\frac{\mathrm{tol}}{e_k}\right)^{1/5}, $$ {#eq-adaptive-step} with safety factor $s<1$, bounded step changes, component scaling and safeguards for zero error. Production solvers combine relative and absolute tolerances; position and angular velocity need not share one unscaled tolerance. Dormand–Prince RK45 is an explicit nonstiff method. [SciPy's solver documentation](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html) distinguishes explicit RK methods from Radau/BDF options for stiff problems. The requested output times are not necessarily the internal step sequence, and event detection can miss multiple crossings within one step. For held inputs, backward Euler and implicit trapezoidal transitions satisfy $$ x_{k+1}=x_k+h f(x_{k+1},u_k,t_{k+1}), $$ {#eq-backward-euler} $$ \begin{aligned} x_{k+1}=x_k+\tfrac h2[& f(x_k,u_k,t_k)\\ &+f(x_{k+1},u_k,t_{k+1})]. \end{aligned} $$ {#eq-trapezoidal} Their nonlinear equations must be solved sufficiently accurately. Differentiate the converged implicit relation, including its inverse Jacobian, when forming transition sensitivities. ### A Convergence Bound With Explicit Assumptions Suppose the trajectory remains in a domain where $$ \|f(x_1,u,t)-f(x_2,u,t)\|\leq L\|x_1-x_2\|, $$ {#eq-lipschitz} and its second time derivative is bounded by $M$. With exact initial state and uniform Euler step $h$, the recursion $e_{k+1}\leq(1+Lh)e_k+Mh^2/2$ yields $$ e_N\leq\frac{M}{2L}(e^{L(T-t_0)}-1)h. $$ {#eq-euler-error-bound} For $L=0$, use the limit $M(T-t_0)h/2$. This is a finite-horizon propagation argument under the stated domain and smoothness assumptions. It does not need a claim that physical perturbations never grow exponentially, or an appeal to a PDE equivalence theorem. For example, Euler converges to the physically unstable solution of $\dot x=x$. Halving $h$ asymptotically halves its fixed-time error. For RK4 the global error ratio approaches 16, whereas the one-step error ratio approaches 32. These are different experiments. ### Absolute Stability and Fast Modes For the scalar test equation $\dot x=\lambda x$, a method has amplification $R(h\lambda)$. Its absolute-stability region satisfies $|R(z)|\leq1$. Euler has $R(z)=1+z$; backward Euler has $R(z)=1/(1-z)$; trapezoidal integration has $R(z)=(1+z/2)/(1-z/2)$. Consider $$ \dot x_1=-0.01x_1,\qquad \dot x_2=-1000x_2. $$ {#eq-stiff-example} Euler requires $h<0.002$ for strict decay of the fast mode. The value $0.002$ is the neutral boundary, not a desirable accuracy choice. Backward Euler and trapezoidal integration are A-stable, but only backward Euler among these two is L-stable. At $z=-100$, their amplifications are $1/101$ and $-49/51$: the latter barely damps an extremely fast physical decay and alternates sign. Implicit does not mean universally stable, and A-stability does not imply accurate resolution of every mode. ### Keep the Error Sources Separate Under suitable bounds, accumulated integration error can be estimated by $$ e_N\leq C N h^{p+1}e^{LNh}, $$ {#eq-total-error} so, for fixed $Nh=T$, $$ e_N=O(h^p). $$ {#eq-fixed-horizon-error} In contrast, for a smooth endpoint map and a finite input perturbation of size $\epsilon$, $$ \mathcal E_{\rm variation}(T,\epsilon)=O(\epsilon^2) $$ {#eq-linearization-residual} when its second derivative is bounded. This error generally approaches a nonzero value as $h\to0$ at fixed $\epsilon$. Part I's $\dot x=x^2$ example demonstrates it exactly. Repeated optimization updates can improve a nominal trajectory, but no universal per-iteration $O(h)$ prediction-error theorem follows. For a swing calculation, report integration convergence, perturbation-amplitude convergence, optimization residuals, constraint/event errors and uncertainty in the mechanical model separately. A tiny integrator tolerance cannot establish that the inferred joint forces or shaft parameters are correct. ### Verify Transition Derivatives The quantities needed by the discrete optimizer are $$ A_k^d=D_xF_k,\qquad B_k^d=D_uF_k. $$ {#eq-jacobian-comp} A central directional check is $$ D_xF_k\,v\approx \frac{F_k(x+\epsilon v,u)-F_k(x-\epsilon v,u)}{2\epsilon}. $$ {#eq-finite-diff-jacobian} Sweep several scales of $\epsilon$ to reveal truncation and roundoff effects, using dimensionally scaled variables. Automatic differentiation differentiates the executed computational graph to floating-point accuracy; it does not remove model error, nondifferentiable branches, solver tolerance effects or a missing event-time derivative. Compare against an independent calculation before interpreting a gradient as physical evidence. ## Control Synthesis via DDP, iLQR, and MPC ### One Quadratic Model for DDP and iLQR Both differential dynamic programming (DDP) and iterative LQR (iLQR) use a quadratic local cost model along a nominal nonlinear rollout. Full DDP also retains second derivatives of the discrete dynamics in the local action-value Hessian; iLQR omits those particular terms. A first-order state-cost expansion is not the defining iLQR approximation. Use $F_k$ for the actual discrete transition and $L_k$ for the discrete running cost, including the chosen quadrature and timestep. For fixed initial state, $$ V_k(x)=\min_{u_k,\ldots,u_{N-1}} \left[\phi(x_N)+\sum_{j=k}^{N-1}L_j(x_j,u_j)\right], $$ {#eq-cost-to-go} subject to $x_{j+1}=F_j(x_j,u_j)$. Bellman's equation is $$ \begin{aligned} V_k(x)&=\min_u[L_k(x,u)+V_{k+1}(F_k(x,u))],\\ V_N&=\phi. \end{aligned} $$ {#eq-bellman} General nonlinear value functions need not be smooth everywhere. The expansions below describe a smooth local model near a chosen nominal trajectory, not a globally quadratic value function. ### Derivatives of the Local Action Value At $\bar x_{k+1}=F_k(\bar x_k,\bar u_k)$, write $$ V_{k+1}(\bar x_{k+1}+\eta) \approx V^+ + (v^+)^\top\eta+\tfrac12\eta^\top P^+\eta. $$ {#eq-quadratic-value} With $A=D_xF_k$, $B=D_uF_k$, $$ \eta=A\delta x+B\delta u+o(\|(\delta x,\delta u)\|). $$ {#eq-thesis-linearized-dynamics} The running-cost expansion retains all quadratic blocks: $$ \begin{aligned} L_k(\bar x+\delta x,\bar u+\delta u)\approx{}& L^0+L_x^\top\delta x+L_u^\top\delta u\\ &+\tfrac12\delta x^\top L_{xx}\delta x +\delta x^\top L_{xu}\delta u\\ &+\tfrac12\delta u^\top L_{uu}\delta u. \end{aligned} $$ {#eq-quadratic-cost} Define the action value using the same absolute-state transition, $$ \mathcal Q_k(x,u)=L_k(x,u)+V_{k+1}(F_k(x,u)). $$ {#eq-q-function} Its local approximation is $$ \begin{aligned} \mathcal Q_k\approx{}&q^0+q_x^\top\delta x+q_u^\top\delta u\\ &+\tfrac12\delta x^\top Q_{xx}\delta x +\delta x^\top Q_{xu}\delta u\\ &+\tfrac12\delta u^\top Q_{uu}\delta u. \end{aligned} $$ {#eq-q-quadratic} The chain rule gives full-DDP derivatives: $$ \begin{aligned} q_x&=L_x+A^\top v^+,\\ q_u&=L_u+B^\top v^+,\\ Q_{xx}&=L_{xx}+A^\top P^+A+\sum_i v_i^+(F_{k,i})_{xx},\\ Q_{xu}&=L_{xu}+A^\top P^+B+\sum_i v_i^+(F_{k,i})_{xu},\\ Q_{uu}&=L_{uu}+B^\top P^+B+\sum_i v_i^+(F_{k,i})_{uu}. \end{aligned} $$ {#eq-q-derivatives} These are derivatives of the discrete transition, not unscaled continuous acceleration derivatives. [Tassa, Mansard and Todorov's DDP treatment](https://roboti.us/lab/papers/TassaICRA14.pdf) describes this distinction between full DDP and its Gauss–Newton variant, as well as the need to treat control limits within the local optimization. ### Solving the Unconstrained Quadratic Subproblem If $Q_{uu}\succ0$ and the input update is unconstrained, $$ q_u+Q_{xu}^\top\delta x+Q_{uu}\delta u=0, $$ {#eq-q-optimality} so $$ \delta u=K\delta x+d, $$ {#eq-optimal-control-update} $$ K=-Q_{uu}^{-1}Q_{xu}^\top,\qquad d=-Q_{uu}^{-1}q_u. $$ {#eq-gains-feedforward} This sign convention also applies to the iLQR update: $$ \delta u_k=K_k\delta x_k+d_k. $$ {#eq-ilqr-control} Use linear solves or factorizations in code, not an explicit matrix inverse. Substitution yields $$ \begin{aligned} V_k^{\rm local}(\delta x) &=q^0-\tfrac12q_u^\top Q_{uu}^{-1}q_u\\ &\quad+v_k^\top\delta x+\tfrac12\delta x^\top P_k\delta x, \end{aligned} $$ {#eq-value-update} $$ \begin{aligned} v_k&=q_x-Q_{xu}Q_{uu}^{-1}q_u,\\ P_k&=Q_{xx}-Q_{xu}Q_{uu}^{-1}Q_{xu}^\top. \end{aligned} $$ {#eq-value-parameters} The local quadratic problem is solved exactly in exact arithmetic under its hypotheses. That does not make its approximation to the nonlinear problem exact or establish a global optimum. ### The iLQR Approximation Dropping only the value-gradient-weighted dynamics Hessians gives $$ \begin{aligned} Q_{xx}&=L_{xx}+A^\top P^+A,\\ Q_{xu}&=L_{xu}+A^\top P^+B,\\ Q_{uu}&=L_{uu}+B^\top P^+B. \end{aligned} $$ {#eq-ilqr-approximation} A scalar check makes the omitted terms visible. Let $F(x,u)=x^2+xu+u^2$, $L(x,u)=\tfrac12x^2+\tfrac12u^2+0.2xu$ and $V^+(y)=\tfrac12y^2$. At $(x,u)=(0.3,-0.2)$, $F=0.07$, $DF=(0.4,-0.1)$ and $D^2F=\begin{pmatrix}2&1\\1&2\end{pmatrix}$. The exact action-value Hessian is $\begin{pmatrix}1.30&0.23\\0.23&1.15\end{pmatrix}$. iLQR omits $0.07D^2F$, while retaining the state and cross curvature of $L$. Independent finite differences verify the full Hessian. ### Complete Iteration and Verification Initialize the nominal states by rolling out chosen inputs from the fixed initial state; initialize $v_N=\phi_x$ and $P_N=\phi_{xx}$. During the backward pass, compute the chosen derivative model and check the input Hessian before factorization. If regularization is needed, state which Hessian or value model is changed and recompute consistently. During the forward pass, use $u_k^{\rm new}=\bar u_k+\alpha d_k+K_k(x_k^{\rm new}-\bar x_k)$ and roll out the nonlinear transition from the original initial state. A line search compares actual cost reduction with the local prediction. Scaling the feedforward step leaves the feedback correction responding to the newly generated states. With active input limits, solve the appropriate constrained quadratic subproblem; arbitrary clipping can fail when the input Hessian couples components. A rollout satisfies its discrete dynamics, but it need not satisfy state constraints, continuous-time limits between nodes or the physical plant. Report stationarity, feasibility, accepted step size and cost reduction. A tiny improvement or an iteration limit can indicate stagnation. ### Convergence Requires More Than Positive Input Curvature Positive $Q_{uu}$ makes a particular local control minimization well posed. It does not alone prove convergence of the nonlinear iteration. Local quadratic Newton convergence requires a nonsingular reduced Hessian, sufficient derivative regularity, a close initialization, and eventually accepted full steps. DDP convergence results depend on the particular algorithm and assumptions. iLQR's Hessian approximation does not generally give superlinear convergence. See the [worked control chapter](../The_Geometry_of_Motion/quarto/ch05_optimal_control.html) for the precise distinction. Use nonlinear rollout, a sufficient-decrease acceptance rule, bounded line search, regularization, and a declared failure exit. If regularization changes the computed policy, evaluate that actual policy in the quadratic model rather than substituting an unregularized Schur-complement formula indiscriminately. Constraints require a compatible constrained subproblem. No benchmark package establishes a universal 5–20 iteration count here. Dense Riccati algebra includes state-control cross products as well as matrix factorizations; derivative and rollout costs depend on the model and implementation. Report timings together with dimensions, horizon, tolerances, derivative method, and achieved residuals. ### Constrained MPC Extension #### Handling Control Bounds For $u_{\min} \leq u \leq u_{\max}$: **Constrained local solve**: Minimize the local quadratic model subject to bounds on the actual updated input. With a positive-definite input Hessian, this is a convex box-constrained quadratic program. An active-set method must account for coupling between free and saturated components. Componentwise clipping is a solution only in special separable cases; the coupled example above shows why. **Quadratic Penalty**: An input-violation penalty can use $\rho \|\max(0, u-u_{\max})\|^2+\rho \|\max(0,u_{\min}-u)\|^2$. This is not an augmented Lagrangian, which also uses multiplier terms and updates. A finite penalty need not enforce exact feasibility, and increasing $\rho$ can worsen conditioning. Check constraint residuals explicitly. Likewise, clipping an unconstrained update is not equivalent to solving the box-constrained quadratic subproblem. #### State Constraints For path constraints $g(x_k) \leq 0$: **Barrier methods**: $$ L_{\text{aug}}(x, u) = L(x, u) - \mu \sum_i \log(-g_i(x)) $$ {#eq-log-barrier} A logarithmic barrier requires a strictly feasible interior. Convergence of a sequence as $\mu\to0$ requires appropriately solved subproblems and regularity assumptions; it is not a general global-optimality guarantee for a nonlinear problem. **Sequential quadratic programming** (SQP): Linearize constraints, solve QP subproblem per iteration. ### Model Predictive Control (MPC) #### Receding Horizon Strategy At each time $t$: 1. **Measure** current state $x(t)$ 2. **Solve** finite-horizon optimal control (via DDP/iLQR): $$\min_{u(\cdot)} \int_t^{t+T} L(x, u) \, d\tau + \phi(x(t+T))$$ {#eq-mpc-problem} 3. **Apply** first control $u^*(t)$ (open-loop for $\Delta t$) 4. **Shift** horizon to $[t+\Delta t, t+T+\Delta t]$, repeat **Receding-horizon replanning** can respond to the following effects, subject to estimation quality, feasibility, and robustness of the chosen design: - Model mismatch - Disturbances - Unmodeled dynamics Replanning uses updated state estimates and a chosen predictive model. Its benefit must be assessed with sensing delay, optimization time, admissibility and model uncertainty included. #### Stability Guarantees ::: {.callout-important} ### MPC Stability Theorem For the sampled model $x^+=F(x,u)$, specify a terminal set $X_f$ and terminal controller $k_f$. Require that $k_f$ respects the constraints, maps $X_f$ into itself, and satisfies $$ \phi(F(x,k_f(x)))-\phi(x)\leq-\ell_d(x,k_f(x)),\qquad x\in X_f, $$ where $\ell_d$ is the stage cost over a sampling interval. With initial feasibility, exact nominal dynamics, and an attained finite-horizon minimum at each step, shifting the previous optimal sequence and appending $k_f$ supplies a feasible next candidate and the value decrease $V_N(x^+)-V_N(x)\leq-\ell_d(x,u^*)$. A suboptimal solver needs a compatible certified-decrease argument; a small optimization residual alone does not prove this inequality. Positive-definite state-cost lower bounds, appropriate local upper bounds and continuity of $V_N$, and admissible bounded trajectories then yield the corresponding stability and attraction result on the certified feasible domain. A nonnegative terminal cost alone is not a Lyapunov certificate, and merely choosing a long horizon does not replace these conditions. Disturbances and model mismatch need a separate robust or practical-stability argument. ::: ### Computational Complexity and Timing Budgets Dense Riccati algebra involves terms of order $N(n^3+n^2m+nm^2+m^3)$, with opportunities to exploit structure. This excludes model-specific derivative evaluation, full-DDP dynamics Hessians, factorization retries, and constrained subproblems. A forward rollout costs $N$ transition evaluations per line-search trial; a general nonlinear transition does not universally cost $O(n)$. Finite differences and automatic differentiation also have different costs depending on the computational graph. At a control rate of 100 Hz, the available sampling interval is 10 ms. That arithmetic is a budget, not evidence that a specified nonlinear optimization meets it. Record hardware, horizon, dimensions, derivative method, worst-case latency, feasibility, and achieved optimization residuals. Warm starts and iteration limits can reduce runtime, but a time limit requires an admissible fallback policy and does not certify convergence. No universal iteration count or real-time benchmark is established here. ### Robust and Stochastic Extensions #### Robust MPC For uncertain dynamics $\dot{x} = f(x, u, w)$ with bounded $\|w\| \leq \epsilon$: **Min-max formulation**: $$ \min_{u(\cdot)} \max_{\|w\| \leq \epsilon} J(x, u, w) $$ {#eq-robust-mpc} **Tube MPC**: Specify a nominal trajectory and a feedback policy whose error set remains bounded under every modeled disturbance; tighten constraints accordingly. An open-loop input sequence and a causal feedback policy define different min-max problems. The uncertainty set, information pattern and invariant-error argument are part of the formulation. #### Stochastic MPC For an Itô model $dx=f(x,u)dt+\sigma(x,u)dW_t$ with Brownian motion $W_t$, admissible controls must be adapted to the available information. One possible risk-sensitive objective is: **Risk-sensitive cost**: $$ \min_{u(\cdot)} \mathbb{E}\left[ \int_0^T L(x, u) \, dt + \phi(x(T)) \right] + \lambda \text{Var}[J] $$ {#eq-risk-sensitive-cost} **Certainty equivalence** has precise validity in standard linear-quadratic-Gaussian settings with their separation assumptions. It is not a general robustness theorem for a nonlinear plant. In general, $\mathbb E[f(x,u)]\ne f(\mathbb E[x],\mathbb E[u])$, and even small noise can matter near constraints, unstable modes, or mode transitions. Evaluate the stochastic or robust formulation actually being used. ### Summary: DDP/iLQR/MPC **Key insights**: 1. **DDP/iLQR** solve local quadratic subproblems whose relation to the nonlinear objective depends on the retained derivatives 2. **Riccati recursion** propagates cost-to-go backward through tangent spaces 3. **Feedback gains** $K_k$ act on local state deviations within the chosen quadratic model 4. **MPC** replans; robustness requires additional design assumptions and verification 5. **Real-time feasibility** depends on measured computation, problem structure, and an admissible fallback These methods use linear algebra and local derivatives in different ways. Their convergence, stability, and performance require the separate assumptions stated above. \newpage ## Case Studies and Worked Examples These examples test different parts of the argument: exact optimization of a stated linear model, configuration-dependent mechanical coupling, and instantaneous underactuation with controllability over time. Reproducible mathematical results are identified explicitly. None is an experimentally validated golf, spacecraft or flight controller. ### Case Study 1: Spacecraft Rendezvous #### Model and Coordinate Convention Use radial $r_x$, along-track $r_y$ and orbit-normal $r_z$ coordinates in the rotating frame of a circular reference orbit. For small relative separation compared with orbital radius, the Clohessy–Wiltshire approximation is $$ \begin{aligned} \ddot r_x-2n\dot r_y-3n^2r_x&=u_x/m,\\ \ddot r_y+2n\dot r_x&=u_y/m,\\ \ddot r_z+n^2r_z&=u_z/m. \end{aligned} $$ {#eq-cw-equations} Here $n$ is the circular reference's orbital rate and $u$ is force in newtons. [MIT's relative-motion derivation](https://ocw.mit.edu/courses/16-346-astrodynamics-fall-2008/e4f0632a9f1c98f7e9b25492e1a30eb1_lec_26.pdf) obtains the linear equations by neglecting higher-order relative-position terms. A constant-coefficient linear model is exact as a mathematical model; its approximation to nonlinear orbital mechanics remains limited by these assumptions. Let $x=(r,\dot r)$, $\dot x=Ax+Bu$, with constant mass. Choose a soft terminal penalty and squared-force cost: $$ J=\tfrac12x(T)^\top Q_fx(T) +\tfrac12\int_0^T u^\top Ru\,dt. $$ {#eq-spacecraft-cost} This objective is not fuel consumption and does not impose exact rendezvous. Weights must include appropriate units or arise from declared nondimensional scales. For the numerical example, evaluate positions in metres, velocities in metres per second, forces in newtons and time in seconds; define the numerical objective with $Q_f=\operatorname{diag}(10^4,10^4,10^4,10^2,10^2,10^2)$ and $R=I$. Thus these coefficients weight numerical quantities in those specified units; changing units requires transforming the weights. #### Reproducible Sampled LQR Result Set $n=0.001$ rad/s, $m=100$ kg, $T=1000$ s and $N=50$ held-input intervals of $h=20$ s. Obtain the exact sampled matrices from $$ \exp\!\left[h\begin{pmatrix}A&B\\0&0\end{pmatrix}\right] =\begin{pmatrix}F&G\\0&I\end{pmatrix}. $$ The discrete running weight is $R_d=hR$. Omitting $h$ changes the objective relative to the terminal penalty. With no running state penalty, the backward recursion is $$ \begin{aligned} P_N&=Q_f,\\ S_k&=R_d+G^\top P_{k+1}G,\\ K_k&=-S_k^{-1}G^\top P_{k+1}F,\\ P_k&=F^\top P_{k+1}F+F^\top P_{k+1}GK_k. \end{aligned} $$ Roll out $u_k=K_kx_k$, $x_{k+1}=Fx_k+Gu_k$ from the actual initial state. This is one Riccati backward pass and one forward rollout, not an iterative nonlinear solver. The following complete calculation uses NumPy and SciPy: ~~~python import numpy as np from scipy.linalg import expm def sampled_cw_example() -> tuple[np.ndarray, np.ndarray, float]: """Solve the stated unconstrained, held-input quadratic problem.""" rate, mass, horizon, intervals = 0.001, 100.0, 1000.0, 50 step = horizon / intervals a = np.zeros((6, 6)) a[:3, 3:] = np.eye(3) a[3, 0], a[3, 4] = 3 * rate**2, 2 * rate a[4, 3], a[5, 2] = -2 * rate, -rate**2 b = np.zeros((6, 3)) b[3:] = np.eye(3) / mass generator = np.zeros((9, 9)) generator[:6, :6], generator[:6, 6:] = a, b transition = expm(step * generator) f, g = transition[:6, :6], transition[:6, 6:] terminal = np.diag([1e4, 1e4, 1e4, 1e2, 1e2, 1e2]) running = step * np.eye(3) p = terminal.copy() gains = [] for _ in range(intervals): gain = -np.linalg.solve(running + g.T @ p @ g, g.T @ p @ f) gains.append(gain) p = f.T @ p @ f + f.T @ p @ g @ gain p = (p + p.T) / 2 states = [np.array([100.0, 50.0, 20.0, 0.01, 0.02, 0.0])] controls = [] for gain in reversed(gains): control = gain @ states[-1] controls.append(control) states.append(f @ states[-1] + g @ control) x, u = np.asarray(states), np.asarray(controls) cost = 0.5 * x[-1] @ terminal @ x[-1] cost += 0.5 * sum(control @ running @ control for control in u) return x, u, float(cost) ~~~ The reproduced terminal position is approximately $(1.66617\times10^{-6},3.04859\times10^{-7},1.48660\times10^{-7})$ m. Its terminal velocity is approximately $(-0.0507127,0.0117398,-0.00839610)$ m/s. The objective is 1.22796775, peak force magnitude 0.09360394 N and accumulated force magnitude 39.39750874 N·s. An independent condensed quadratic-program calculation agrees with the Riccati control sequence to approximately $10^{-8}$ N; exact sampled propagation also agrees with independently integrated continuous dynamics. The very small position error coexists with substantial terminal velocity. That is the actual tradeoff of these soft weights, not successful docking. Exact endpoint constraints, thrust limits, collision avoidance, actuator allocation and nonlinear orbital-model verification would define a different problem. For an ideal thruster at fixed specific impulse $I_{\rm sp}$, propellant relates to delivered impulse through $m_p=I_{\rm total}/(g_0I_{\rm sp})$, as follows from [NASA's specific-impulse definition](https://www1.grc.nasa.gov/beginners-guide-to-aeronautics/specific-impulse/). For example, 12.4 N·s at 300 s corresponds to **0.00421483 kg**, using $g_0=9.80665$ m/s². This separate unit check is not the impulse of the simulation above. Real thruster allocation can make the sum of individual thrust magnitudes exceed the norm of the resultant, and squared force is a different integral again. ### Case Study 2: Planar Robot Arm #### Placement Determines Both Gravity and Kinematics Choose two massless links of length 1 m, with a 1 kg point mass at each distal end. This is not the inertia model of two uniform massive rods. Let $q_1$ be measured counterclockwise from the horizontal world $+x$ axis and $q_2$ be the relative elbow angle. Then $$ \begin{aligned} p_1&=(\cos q_1,\sin q_1),\\ p_2&=p_1+(\cos(q_1+q_2),\sin(q_1+q_2)). \end{aligned} $$ Writing $$ M(q)\ddot q+C(q,\dot q)\dot q+G(q)=\tau, $$ {#eq-thesis-robot-dynamics} the chosen masses and lengths give $$ M(q)= \begin{pmatrix} 3+2\cos q_2&1+\cos q_2\\ 1+\cos q_2&1 \end{pmatrix}, $$ {#eq-inertia-matrix} in kg·m². With $a=-\sin q_2$ in the corresponding mass-length units, $$ C(q,\dot q)= \begin{pmatrix} a\dot q_2&a(\dot q_1+\dot q_2)\\ -a\dot q_1&0 \end{pmatrix}. $$ {#eq-coriolis-matrix} The potential is $U=g[2\sin q_1+\sin(q_1+q_2)]$ in joules with these SI parameter values, so $$ G(q)=g \begin{pmatrix} 2\cos q_1+\cos(q_1+q_2)\\ \cos(q_1+q_2) \end{pmatrix}. $$ {#eq-gravity-vector} Gravity contains cosines for angles measured from the horizontal. A sine formula would require a different angular origin and matching kinematics. Independent checks construct $M$ from the two point-mass Jacobians, differentiate $U$ numerically and verify $\dot q^\top(\dot M-2C)\dot q=0$. Thus the mechanical energy satisfies $\dot E=\dot q^\top\tau$. These checks link placement, inertia, gravity and power; checking only that $M$ is positive definite would miss the inconsistent gravity convention. #### Inverse Kinematics Before Trajectory Optimization For a target $(x,y)$ with general positive link lengths, define $c_2=(x^2+y^2-L_1^2-L_2^2)/(2L_1L_2)$. Reachability requires $|L_1-L_2|\leq\sqrt{x^2+y^2}\leq L_1+L_2$. Reject targets outside this annulus apart from a stated numerical tolerance; clipping an arbitrary unreachable target into $[-1,1]$ invents a different target. For the branch $q_2=+\arccos c_2$, take $q_1=\operatorname{atan2}(y,x)-\operatorname{atan2}(L_2\sin q_2,L_1+L_2\cos q_2)$. The negative branch gives the other elbow configuration when distinct. At a singular folded configuration, additional care is needed because the shoulder angle can be indeterminate. With unit lengths, the checked angles are: - Target $(0.5,0.5)$ m: $q\approx(-24.29518895^\circ,138.59037789^\circ)$. - Target $(0.3,0.8)$ m: $q\approx(4.73396255^\circ,129.41998447^\circ)$. By contrast, $(30^\circ,60^\circ)$ places the end at $(\sqrt3/2,1.5)$ m. Forward substitution verifies each result. These are kinematic endpoints, not a computed dynamically feasible trajectory or a convergence benchmark. #### What Coupling Establishes The off-diagonal entries of $M$ mean that a torque at one joint can accelerate both joints. Their dependence on posture changes that instantaneous mapping. It does not mean the tangent-space dimension changes, and matrix entries do not define a universal ordering of configurations from least to greatest inertia. DDP can use derivatives along a changing nominal motion, while a fixed linear controller can still work within a verified operating region. Neither possibility is decided solely by observing that $M(q)$ varies. For a golf model, anatomical constraints, shaft modes, ground reactions and available actuator directions add further structure; the fully actuated two-link arm does not establish a universal optimal sequencing rule. ### Case Study 3: Quadrotor Stabilization #### Hover in a Consistent Frame Let world and body $e_3$ point upward at level attitude, and let $R$ map body vectors to world vectors. Rotor thrusts $u_i\geq0$ act along body $+e_3$. Let $q$ be a unit, scalar-first quaternion for $R$ and $\omega$ be body angular velocity. Then $$ \begin{aligned} \dot p&=v,\\ \dot v&=-g e_3+\frac{\sum_i u_i}{m}R e_3,\\ \dot q&=\tfrac12q\otimes(0,\omega),\\ \dot\omega&=J^{-1}(\tau-\omega\times J\omega). \end{aligned} $$ {#eq-quadrotor-dynamics} The pure quaternion $(0,\omega)$ is essential to the quaternion product. This continuous equation preserves $\|q\|=1$; a numerical method must also control its constraint error. Place rotors at body locations $(L,0,0)$, $(0,L,0)$, $(-L,0,0)$ and $(0,-L,0)$. A consistent alternating-spin allocation is $$ \begin{aligned} \tau_x&=L(u_2-u_4),\\ \tau_y&=L(u_3-u_1),\\ \tau_z&=c(u_1-u_2+u_3-u_4). \end{aligned} $$ Here $u_i$ is thrust in newtons, $L$ is in metres and $c$ is a torque-to-thrust ratio in metres. These inputs are not raw motor speed commands: if $u_i=k_f\Omega_i^2$, the model is affine in thrust but nonlinear in $\Omega_i$. At $R=I$, zero rates and any fixed position, $u_i=mg/4$ balances gravity and produces zero moment. Adding gravity and upward thrust with the same sign would not admit this hover equilibrium. #### Twelve Local Errors for a Thirteen-Component State The embedded vector $(p,v,q,\omega)$ has 13 entries but lies on a 12-dimensional state manifold. Near a fixed hover reference, use the stacked error $$ e=\begin{pmatrix} p-p_{\rm ref}\\ v\\ \operatorname{Log}(R_{\rm ref}^\top R)^\vee\\ \omega \end{pmatrix}. $$ The logarithm is a local rotation coordinate, with branch restrictions away from the identity; it is not a global nonsingular three-parameter attitude representation. A $12\times12$ quadratic weight acts on this 12-component error, not on the raw 13-component subtraction. The construction also respects that $q$ and $-q$ represent the same attitude. A running cost can penalize $e$ and $\delta u=u-u_{\rm hover}$ so the desired steady operation has zero deviation cost. Moving references require consistent angular-rate transport and tracking errors. #### Instantaneous Rank and Controllability Over Time For $m=1.2$ kg, $L=0.2$ m, $c=0.01$ m and $J=\operatorname{diag}(0.02,0.02,0.04)$ kg·m², the linearized input matrix has rank four. At level hover, $\delta\dot v_x=g\,\delta\theta_y$ and $\delta\dot v_y=-g\,\delta\theta_x$: rotor-induced moments first change attitude, and attitude then changes horizontal acceleration. The controllability matrix $[B,AB,\ldots,A^{11}B]$ has rank 12 for this example. This does not create twelve independent instantaneous inputs. It shows that the linearized state can be controlled over time near an interior hover input, subject to the local-model assumptions. Positive thrust bounds and available thrust margins still constrain finite maneuvers. At 40 Hz, 50 intervals span 1.25 s; a two-second horizon requires 80 intervals. These are sampling facts. Settling time, centimetre tracking accuracy and robustness percentages require an implemented controller, declared perturbations and actual verification; no such performance claim is established here. ### What the Examples Contribute to Swing Analysis The orbital example shows why an exactly solved quadratic objective can still produce an unacceptable terminal event. The arm shows that kinematics, inertia, gravity and energy must share one physical placement. The quadrotor shows why instantaneous input rank and the ability to influence a future task are different questions. Together they motivate a disciplined swing argument: derive a consistent mechanical model, specify the available inputs and intervention, propagate the resulting motion, and evaluate the actual contact objective with its constraints and uncertainty. Local derivatives make parts of this process tractable; they do not replace the physical and experimental checks. ### Reproduction Scope The companion [control and case-study checks](https://github.com/D-sorganization/AffineDrift/blob/main/tests/test_tangent_control_rigor.py) reproduce the Hamiltonian signs, constrained quadratic example, integration orders, DDP curvature, sampled orbital calculation, arm mechanics and hover rank calculations. They are verification examples rather than a released nonlinear-control package or empirical validation. ## Conclusion and Future Directions ### Three Levels of Mechanical Reasoning The useful result is a hierarchy of statements with different scopes. 1. **An immediate acceleration change.** In a model affine in the declared force input, the increment is $M^{-1}\mathcal B\,\Delta u$. State, time, model, parameters and constraint mode remain fixed. 2. **A small change in a later output.** The derivative of the solution/output map propagates through the variational equation. Reference motion, perturbation definition, policy and endpoint convention remain fixed. 3. **A finite intervention.** A new nonlinear trajectory determines its own task and event outputs. The intervention and comparison protocol are specified, while the state evolves independently. A tangent space supplies a place for velocity vectors. Differentiation supplies a linear map between tangent spaces. Neither fact makes the underlying vector field or its finite flow linear. First variations superpose because they obey one linear equation along a common reference, not because all nonlinear motion is a sequence of exactly linear physical systems. ### How the Swing's Components Interact At one posture, an applied joint torque can accelerate several coordinates through the inverse inertia matrix. Shaft deflection changes the geometry and energy stores entering that matrix and the other loads. Subsequent posture and velocity changes alter both the input gain and the reference transport. Finally, clubhead kinematics and contact timing select which of those state changes matter for the strike. This sequence connects local mechanical coupling to an eventual task outcome without identifying them. A large instantaneous acceleration contribution can have a small endpoint effect. Positive input work can accompany a decrease in a particular speed. A small interaction residual between two perturbations can coexist with substantial model bias or a large single-trajectory approximation error. Input definition also matters. Replaying a torque history, preserving a feedback policy, changing stiffness, changing excitation in an activation model, and reoptimizing after an equipment change are different experiments. A result from one should not be described as the result of another. ### What Linear Control Can Establish For an equilibrium with the required smoothness, a Hurwitz closed-loop Jacobian can support a local nonlinear stability theorem. LQR additionally solves a specified linear-quadratic optimization problem under its stabilizability and detectability conditions. Finite-horizon trajectory control, DDP/iLQR and MPC need their own conditions on feasible references, objectives, constraints, numerical rollout and closed-loop execution. An exact derivative is valuable input to those methods. It does not by itself establish global optimality, a region of attraction, a guaranteed contact result or a biological explanation. A certificate must state the model, domain, disturbance assumptions and physical output it controls. ### A Verification Protocol for the Argument For a proposed golf application, first specify the state, coordinate convention, actuator inputs, constraints and measured outputs. Demonstrate dimensional consistency and force/energy balances. Then test derivative code against independent finite perturbations on smooth intervals. Refine the integration tolerance and perturbation amplitude separately. Compare outcomes at the same declared time or event. Include event-time and reset sensitivity when appropriate. Test the actual proposed feedback or input replay, retain correlated uncertainty, and report the output scale and tolerance used to judge a discrepancy. Verify a finite nonlinear intervention directly before presenting its first-order estimate as quantitatively adequate. Finally, distinguish numerical verification from empirical validation. An internally consistent example can test an equation or algorithm. Evidence for a golfer requires identified parameters, measured loads and motion, uncertainty analysis, and prediction on data not used to fit the model. Local mathematical correctness is necessary for that investigation, but cannot substitute for it. ### Open Technical Questions The next research steps are concrete: identify which input directions can be estimated from available measurements; determine the domain over which finite predictions meet task tolerances; and test how contact transitions, flexible modes, actuator limits and feedback alter those conclusions. A proposed geometric residual bound must distinguish coordinate Hessians, dynamic transport and intrinsic curvature. A stochastic extension needs a specified stochastic differential equation and calculus convention; Brownian paths do not have ordinary instantaneous velocities to which the deterministic argument can simply be applied. These questions do not have universal answers based on swing duration or the dimension of a tangent space. They can be investigated through explicit models, controlled comparisons and evidence. ::: {.callout-note} ### In Plain Language A snapshot can tell us how an extra torque changes acceleration at the current posture. Predicting the rest of the swing requires following how that change alters the body's motion, the shaft and the later effect of every load. A small-change calculation is a precise way to study that connection near a reference swing. A larger change needs a new simulation and, for claims about real golfers, measurements that test its prediction. ::: ## References {.unnumbered} 1. Abraham, R., & Marsden, J. E. (1978). *Foundations of Mechanics* (2nd ed.). Benjamin/Cummings. 2. Arnold, V. I. (1989). *Mathematical Methods of Classical Mechanics* (2nd ed.). Springer. 3. Bryson, A. E., & Ho, Y.-C. (1975). *Applied Optimal Control*. Hemisphere Publishing. 4. Jacobson, D. H., & Mayne, D. Q. (1970). *Differential Dynamic Programming*. Elsevier. 5. Khalil, H. K. (2002). *Nonlinear Systems* (3rd ed.). Prentice Hall. 6. Lee, J. M. (2012). *Introduction to Smooth Manifolds* (2nd ed.). Springer. 7. Liberzon, D. (2003). *Switching in Systems and Control*. Birkhäuser. 8. Mayne, D. Q., et al. (2000). Constrained model predictive control: Stability and optimality. *Automatica*, 36(6), 789-814. 9. Sastry, S. (1999). *Nonlinear Systems: Analysis, Stability, and Control*. Springer. 10. Tassa, Y., et al. (2012). Synthesis and stabilization of complex behaviors through online trajectory optimization. In *Proc. IEEE/RSJ IROS* (pp. 4906-4913). \newpage ## Appendices {.unnumbered} ## Appendix A: Mathematical Prerequisites A smooth manifold is described locally by coordinates with smooth invertible transition maps. $T_x\mathcal M$ contains curve velocities at one point; $T\mathcal M$ collects these spaces over the manifold. A covector is a linear functional on a tangent space. Generalized force and the differential of a scalar cost are covectors, whereas velocity and a first state variation are vectors. A smooth map has a derivative between the appropriate tangent spaces. Coordinates represent this derivative by a Jacobian. A metric supplies lengths and identifies vectors with covectors, while a connection supplies a rule for differentiating vector fields along curves. Neither structure is required merely to define a tangent space or the derivative of a smooth map. If used, each must be specified. For rotations, use a local chart or a declared Lie-group representation with its frame convention. A coordinate singularity is not automatically a physical instability. For a mechanical state, distinguish configuration from its velocity state, and distinguish both from a clubhead output. These definitions are developed with physical examples in the [configuration-manifold chapter](../The_Geometry_of_Motion/quarto/volume2.qmd#configuration-manifolds). ## Appendix B: Proofs and Derivations This appendix supplies explicit bounds and convergence statements used in the main text. ## B.1 A Finite-Difference Bound From the Endpoint Map Let $F(a,b)=x(t_1;\bar u+a v_1+b v_2)$ use the same initial state for every pair of amplitudes. Assume a $C^2$ solution map exists on a rectangle containing $0\leq a,b\leq\epsilon$. The fundamental theorem of calculus, applied twice componentwise, gives the exact identity $$ r=F(\epsilon,\epsilon)-F(\epsilon,0)-F(0,\epsilon)+F(0,0) =\int_0^\epsilon\int_0^\epsilon F_{ab}(a,b)\,db\,da. $$ Consequently, $$ \|r\|\leq\epsilon^2\sup_{a,b\in[0,\epsilon]}\|F_{ab}(a,b)\|. $$ The mixed second-variation equation in the main text computes $F_{ab}$, using derivatives along each corresponding trajectory. A bound on the state Hessian alone omits mixed state-input terms and nonlinear input effects. For control-affine dynamics $f_0(x)+G(x)u$, the term $D G(x)[\xi]v$ remains even though the pure input Hessian vanishes. Existence of one nominal trajectory is also insufficient: all perturbations used by the bound must remain in the smooth model domain. ## B.2 Intrinsic Curvature Is a Different Quantity For a connection $\nabla$, the curvature identity is $$ R(v,w)\xi=\nabla_v\nabla_w\xi-\nabla_w\nabla_v\xi-\nabla_{[v,w]}\xi. $$ The final Lie-bracket term vanishes for commuting coordinate vector fields, but cannot generally be omitted. Curvature can govern parallel-transport holonomy. That geometric operation is not the same as comparing four forced trajectories of an arbitrary vector field. The counterexample $\dot x=u^2$ lives on the flat Euclidean line, where $R=0$, yet its input interaction is nonzero. It disproves a general bound of this interaction solely by sectional curvature. Coordinate-dependent field derivatives, a specified connection, and intrinsic manifold curvature must be distinguished before proposing a geometric interpretation. ## B.3 Integration Error and Linearization Error For explicit Euler applied to a sufficiently smooth initial-value problem over a fixed finite interval, bounded local truncation error $O(h^2)$ and a suitable Lipschitz stability bound yield global integration error $O(h)$. The constant depends on the interval and growth of perturbations. Lipschitz continuity alone does not supply the higher differentiability needed for every claimed high-order method. If the integrated model is the first-variation approximation rather than the full nonlinear system, reducing $h$ converges to that approximate model. It does not remove finite-perturbation error. For a single smooth input direction, write $x_\epsilon=\bar x+\epsilon\xi+\tfrac12\epsilon^2\eta+o(\epsilon^2)$. A computed first-variation endpoint with a uniform Euler error bound satisfies $$ \|x_\epsilon(t_1)-x_{\mathrm{linear},h}(t_1)\| \leq C_{\mathrm{disc}}h+\tfrac12\epsilon^2\|\eta(t_1)\|+o(\epsilon^2). $$ The self-variation $\eta$ includes state, input, and mixed derivatives and their transport. This distinguishes convergence of the numerical solver from fidelity of the linearized physical prediction. Assess the two with separate step-size and perturbation-amplitude refinements. ## Appendix C: Notation Guide and Glossary ## Primary Notation | Symbol | Meaning | |--------|---------| | $x \in \mathbb{R}^n$ | State vector | | $u \in \mathbb{R}^m$ | Control input | | $f: \mathbb{R}^n \times \mathbb{R}^m \to \mathbb{R}^n$ | Vector field (dynamics) | | $A = \frac{\partial f}{\partial x}$ | State Jacobian | | $B = \frac{\partial f}{\partial u}$ | Input Jacobian | | $\delta x$, $\delta u$ | Infinitesimal perturbations | | $\Phi(t_1, t_0)$ | State transition operator | | $\lambda(t)$ | Costate (adjoint variable) | | $H(x,u,\lambda,t)$ | Hamiltonian | | $L(x,u,t)$ | Optimization running cost; distinct from a mechanical Lagrangian | | $r$ | Residual (failure of superposition) | ## Glossary **Costate** ($\lambda$): A covector propagated by an adjoint equation for a specified cost and comparison protocol. An open-loop future-cost sensitivity need not equal the gradient of an optimized value function. **Fréchet Derivative**: A bounded linear map whose remainder is little-o of the perturbation norm. Existence of derivatives along separate lines is insufficient. **Jacobian**: Coordinate matrix of a differentiable map's derivative. A state-rate Jacobian along a moving coordinate transformation includes the associated basis-change term. **Manifold**: Space that locally resembles Euclidean space; may be curved globally. **Residual**: Deviation from additivity when applying multiple perturbations; $O(\epsilon^2)$ interaction under stated smoothness and fixed-horizon assumptions. **State Transition Operator** ($\Phi$): Linear map transporting perturbations forward in time via variational dynamics. **Tangent Space** ($T_x \mathcal{M}$): Linear space of velocity vectors at point $x$ on manifold $\mathcal{M}$. **Variational Equation**: The exact evolution equation for the derivative of a trajectory family with respect to perturbation amplitude, along one common reference. Its finite-amplitude use is an approximation. ## Appendix D: Code Repository Index The companion regression file [test_tangent_foundation_rigor.py](https://github.com/D-sorganization/AffineDrift/blob/main/tests/test_tangent_foundation_rigor.py) independently checks the differentiability counterexamples, coordinate transformation, pendulum energy, ordered transport, pulse kernel, functional/adjoint sensitivity and combined-motion work example. The existing [residual checks](https://github.com/D-sorganization/AffineDrift/blob/main/tests/test_tangent_residual_rigor.py) check the evolved baseline and mixed-variation transport. These are verification examples, not a released general-purpose nonlinear-control package. The [control and case-study checks](https://github.com/D-sorganization/AffineDrift/blob/main/tests/test_tangent_control_rigor.py) independently verify Part III's optimization, discretization, orbital, arm and hover examples. The sampled orbital calculation is reproducible; no empirical flight or golf performance is claimed. --- **End of Thesis**

[1] Abraham, R., & Marsden, J. E. (1978). Foundations of Mechanics (2nd ed.). Benjamin/Cummings.

[2] Arnold, V. I. (1989). Mathematical Methods of Classical Mechanics (2nd ed.). Springer.

[3] Bryson, A. E., & Ho, Y.-C. (1975). Applied Optimal Control. Hemisphere Publishing.

[4] Jacobson, D. H., & Mayne, D. Q. (1970). Differential Dynamic Programming. Elsevier.

[5] Khalil, H. K. (2002). Nonlinear Systems (3rd ed.). Prentice Hall.

[6] Lee, J. M. (2012). Introduction to Smooth Manifolds (2nd ed.). Springer.

[7] Liberzon, D. (2003). Switching in Systems and Control. Birkhäuser.

[8] Mayne, D. Q., et al. (2000). Constrained model predictive control: Stability and optimality. Automatica, 36(6), 789-814.

[9] Sastry, S. (1999). Nonlinear Systems: Analysis, Stability, and Control. Springer.

[10] Tassa, Y., et al. (2012). Synthesis and stabilization of complex behaviors through online trajectory optimization. In Proc. IEEE/RSJ IROS (pp. 4906-4913).

Part I: Geometric Foundations and Local Dynamics

1 Introduction and Motivation

1.1 The Central Thesis

1.2 Central Claim

1.3 Scope: No Trajectory Superposition

1.4 Terminology: “Tangent Hyperplane” as Affine Subspace

1.5 Scope and Prerequisites

2 Tangent Spaces and Exact Differentiation

2.1 Nonlinear Dynamics and Smoothness

2.2 Continuous Differentiability Does Not Give a Quadratic Bound

2.3 Directional Derivatives

2.4 The Jacobian Matrices

2.5 A Tangent Space Is Not a Derivative

2.6 Coordinate Changes Along a Moving Motion

2.7 Example: The Simple Pendulum

2.8 Exact Superposition of Input Increments

2.9 Counterfactual Decomposition

3 Time Evolution and Moving Tangent Spaces

3.1 Why Global Superposition Fails

3.2 Flow Transport Is Not Parallel Transport

3.3 Separate the Two Refinement Experiments

3.4 Variational Equation

4 Local Optimal Control: LQR and Lyapunov Stability

4.1 Connection to Linear Quadratic Regulation

4.1.1 The LQR Problem

4.1.2 What a Nonlinear LQR Design Establishes

4.1.3 Gain Scheduling Requires Its Own Check

4.2 Lyapunov Stability Analysis

4.2.1 Lyapunov’s Indirect Method

4.2.2 Why This Works

4.2.3 Design and Validation

5 Residuals, Reference Motion, and Nonlinearity

5.1 An Interaction Requires Four Comparable Trajectories

5.2 First Variations and the Mixed Second Variation

5.3 Transport Is Part of the Bound

5.4 What a Small Residual Does Not Establish

Part II: Integral Superposition and Transport

6 Fundamentals of Integral Accumulation

6.1 Integration Preserves Linearity

6.2 The Integration Theorem

6.3 What Superposes: Variations About One Reference

7 State Transition Operator

7.1 Solution of the Variational Equation

7.2 Interpretation: Transport Then Accumulate

7.3 Time-Invariant Special Case

8 Variational Dynamics and Sensitivity Analysis

8.1 Sensitivity to Initial Conditions

8.2 An Input Kernel Is a Density

8.3 The Endpoint Derivative Is an Operator

9 Adjoint (Costate) Dynamics

9.1 Backward Sensitivity

9.2 Adjoint Evolution Equation

9.3 The Policy Determines the Sensitivity

9.4 A Reproducible Gradient Check

9.5 Task and Event Interpretation

10 Physical Interpretations: Impulse and Work

11 Force-Momentum: Impulse Integral

11.1 Newton’s Second Law

11.2 Impulse Definition

11.3 Superposition of Impulses

12 Power-Energy: Work Integral

12.1 Power and Work

12.2 Superposition of Work

13 General Principle: Integral Accumulation

14 Finite-Horizon Residuals and Error Quantification

14.1 Keep the Baseline and the Perturbation Definition Fixed

14.2 Interaction and Single-Trajectory Approximation Error

14.3 A Task-Based Interpretation for the Swing

Part III: Hamiltonian Structure and Optimal Control

15 Variational Principles and Hamiltonian Formulation

15.1 Mechanics and Optimization Ask Different Questions

15.2 A Consistent Minimization Convention

15.3 Costate as Sensitivity

15.4 Deriving the Dynamic Residual Correctly

15.5 Endpoint and Time Conditions

15.6 Bounded Inputs and Bang-Bang Control

15.7 Minimum Squared Input Is Not Mechanical Energy

15.8 A Complete Minimum-Time Example

15.9 LQR as a Special Boundary-Value Problem

15.10 Box Constraints Require the Correct Metric