Trajectory Sheaves

Module for trajectory-valued sheaves encoding linear dynamical system state evolution, together with colocation boundary-value solvers via harmonic extension.

This module provides two trajectory-sheaf constructions:

Autonomous (TrajectorySheaf)

A TrajectorySheaf wraps an EuclideanSheaf on a path graph of length k+1 whose restriction maps encode the autonomous discrete-time dynamics $x_{t+1} = A x_t$. Global sections of the inner sheaf correspond exactly to valid k-step trajectories.

• TrajectorySheaf — struct holding k, A, and the inner sheaf.
• build_trajectory_sheaf — construct the path-graph sheaf for given F, A, and k.
• colocation_trajectory — solve the two-point boundary-value problem $x_0 = x_{\text{given}}$, $x_k = x_{\text{given}}$ via harmonic extension.

Controlled (ControlledTrajectorySheaf)

A ControlledTrajectorySheaf wraps an EuclideanSheaf on a path graph of length k+2 whose restriction maps encode the zero-order-hold discretization of the continuous-time LTI system $\dot{x}(t) = A_c x(t) + B_c u(t)$. Global sections of the inner sheaf correspond exactly to feasible sampled state-control trajectories.

• continuous_to_discrete_zoh — discretize a continuous-time LTI system using the zero-order-hold block-exponential formula.
• ControlledTrajectorySheaf — struct and constructor; builds the inner path-graph sheaf eagerly.
• feasible_control_trajectory_basis — return a particular feasible trajectory and a null-space basis for endpoint-preserving perturbations.

ControlledTrajectorySheaf{T}

Trajectory sheaf for the zero-order-hold discretization of the continuous-time controlled system ẋ(t) = Ac*x(t) + Bc*u(t).

The inner sheaf is defined on a path graph with k+2 vertices:

• vertex 1 is a dummy initial vertex storing only x₁,
• vertices 2, …, k+1 store (xₜ, uₜ) for t = 1, …, k,
• vertex k+2 stores only x_{k+1}.

Global sections of sheaf correspond exactly to feasible sampled state-control trajectories of the discretized system.

Fields:

• k :: Int — number of time steps.
• h :: T — sample period.
• Ac :: AbstractMatrix{T} — continuous-time state matrix.
• Bc :: AbstractMatrix{T} — continuous-time input matrix.
• Ad :: Matrix{T} — discrete-time state matrix.
• Bd :: Matrix{T} — discrete-time input matrix.
• sheaf :: EuclideanSheaf{T} — inner path-graph sheaf.
• state_dim :: Int — state dimension n.
• control_dim :: Int — control dimension m.

ControlledTrajectorySheaf(F::EuclideanSheaf{T},
                          Ac::AbstractMatrix{T},
                          Bc::AbstractMatrix{T},
                          h::Real,
                          k::Int) where T
    -> ControlledTrajectorySheaf{T}

Construct a ControlledTrajectorySheaf from a base sheaf F, continuous-time system matrices (Ac, Bc), sample period h, and number of steps k.

The state dimension n is the total 0-cochain dimension of F (i.e. sum(vertex_stalks(F))). The inner sheaf has k+2 vertices:

• vertex 1: stalk $\mathbb{R}^n$ (dummy initial state),
• vertices $2, \ldots, k+1$: stalk $\mathbb{R}^{n+m}$ (state and control),
• vertex $k+2$: stalk $\mathbb{R}^n$ (terminal state).

The $k+1$ edges enforce the zero-order-hold discrete dynamics:

• dummy initial edge $(1,2)$: $x_1^{\mathrm{dummy}} = x_1^{\mathrm{traj}}$,
• dynamics edges $(t+1, t+2)$ for $t = 1, \ldots, k-1$: $A_d x_t + B_d u_t = x_{t+1}$,
• final dynamics edge $(k+1, k+2)$: $A_d x_k + B_d u_k = x_{k+1}$.

Throws

• ArgumentError if k < 1.
• ArgumentError if h <= 0.
• ArgumentError if Ac is not square.
• ArgumentError if size(Bc, 1) ≠ sum(vertex_stalks(F)).

TrajectorySheaf{T}

A trajectory-valued sheaf encoding k steps of the discrete-time linear dynamics $x_{t+1} = A x_t$.

Fields:

• k :: Int — number of time steps.
• A :: AbstractMatrix{T} — state-transition matrix ($d \times d$).
• sheaf :: EuclideanSheaf{T} — inner path-graph sheaf; its global sections are exactly the valid trajectories $[x_0, x_1, \ldots, x_k]$.

TrajectorySheaf(F::EuclideanSheaf{T}, A::AbstractMatrix{T}, k::Int)
    -> TrajectorySheaf{T}

Construct a TrajectorySheaf from a base sheaf F, dynamics matrix A, and number of steps k.

The inner sheaf is built eagerly by build_trajectory_sheaf.

Throws

• ArgumentError if k < 1.
• ArgumentError if A is not square or has the wrong size relative to F.

nullspace_trajectory_family(
    ts::ControlledTrajectorySheaf{T},
    z_p::AbstractVector{T},
    null_basis::AbstractMatrix{T};
    amplitude::Real = one(T),
    include_negative::Bool = false
) where T
    -> Vector{BlockVector{T}}

Construct a family of feasible trajectories by activating one nullspace basis direction at a time in the affine parameterization z = z_p + null_basis * α.

For each basis column n_j, this returns z_p + amplitude*n_j. If include_negative=true, it also returns z_p - amplitude*n_j.

The returned trajectories use the public block layout: k+1 state blocks of size ts.state_dim followed by k control blocks of size ts.control_dim.

build_trajectory_sheaf(F::EuclideanSheaf{T}, A::AbstractMatrix{T}, k::Int)
    -> EuclideanSheaf{T}

Construct the inner path-graph sheaf for a TrajectorySheaf.

The state dimension d equals the total dimension of the 0-cochain space of F (i.e. sum(vertex_stalks(F))). The returned sheaf has k+1 vertices, each with stalk $\mathbb{R}^d$, connected by k edges with an induced orientation in the coboundary map. For each edge between vertices t and t+1, using the induced orientation $t \to t+1$:

• The restriction map from vertex t is $A$ (dynamics matrix).
• The restriction map from vertex t+1 is $I_d$ (identity).

A 0-cochain $[x_0, \ldots, x_k]$ is a global section iff $A x_t - x_{t+1} = 0$, i.e. $x_{t+1} = A x_t$ (forward dynamics),

Arguments

• F — base EuclideanSheaf whose total 0-cochain dimension gives d.
• A — $d \times d$ state-transition matrix.
• k — number of time steps (must be ≥ 1).

Throws

• ArgumentError if k < 1.
• ArgumentError if A is not square.
• ArgumentError if size(A, 1) ≠ d.

colocation_trajectory(tsheaf::TrajectorySheaf{T}, x0, xk)
    -> BlockVector

Solve the two-point boundary-value problem for the trajectory sheaf by harmonic extension.

Given initial state x0 (assigned to time step 1) and final state xk (assigned to time step k+1), returns the minimum-energy 0-cochain of tsheaf.sheaf with those boundary values. If the boundary data are consistent with the dynamics over k steps (that is, xk ≈ A^k * x0), this coincides with the exact trajectory determined by those endpoints. If the boundary data are inconsistent, there is no exact trajectory with those endpoints; in that case the returned value is the harmonic minimum-energy / least-squares 0-cochain and generally will not satisfy x[t+1] = A * x[t] exactly at every step.

The returned BlockVector has k+1 blocks where trajectory[Block(t)] gives the state at time step t-1 (so Block(1) is x_0, Block(2) is x_1, …, Block(k+1) is x_k).

Arguments

• tsheaf — a TrajectorySheaf{T}.
• x0 — initial state vector of length d = size(tsheaf.A, 1).
• xk — final state vector of length d.

Throws

• ArgumentError if length(x0) ≠ d or length(xk) ≠ d.

continuous_to_discrete_zoh(Ac::AbstractMatrix{T},
                           Bc::AbstractMatrix{T},
                           h::Real) where T
    -> (Ad::Matrix{T}, Bd::Matrix{T})

Discretize the continuous-time LTI system ẋ(t) = Ac*x(t) + Bc*u(t) using zero-order hold on a uniform step size h.

The returned matrices satisfy x[k+1] = Ad*x[k] + Bd*u[k].

The computation uses the block-matrix exponential identity

\[\exp\!\left(h \begin{bmatrix} A_c & B_c \\ 0 & 0 \end{bmatrix}\right) = \begin{bmatrix} A_d & B_d \\ 0 & I_m \end{bmatrix}.\]

Arguments

• Ac — $n \times n$ continuous-time state matrix.
• Bc — $n \times m$ continuous-time input matrix.
• h — sample period (must be positive).

Throws

• ArgumentError if h <= 0.
• ArgumentError if Ac is not square.
• ArgumentError if size(Bc, 1) ≠ size(Ac, 1).

controlled_vehicle_platoon_with_herding(
    Ac::AbstractMatrix{T},
    Bc::AbstractMatrix{T},
    h::Real,
    k::Int,
    herding_gain::Real;
    k_herding::Real = herding_gain
) -> ControlledTrajectorySheaf{T}

Construct a ControlledTrajectorySheaf for a vehicle platoon where herding behavior emerges from the sheaf Laplacian rather than being explicitly encoded in the dynamics.

The resulting sheaf is defined on a ladder graph with:

• Two parallel paths (backbones) representing each vehicle's state trajectory over time
• Each backbone has k+2 vertices following the ControlledTrajectorySheaf pattern
• Rungs connect corresponding time steps between vehicles
• Herding gain controls the strength of attraction between vehicles

Arguments

• Ac: n×n continuous-time state matrix for each vehicle's individual dynamics
• Bc: n×m continuous-time input matrix for each vehicle
• h: sample period for discretization
• k: number of time steps in the trajectory
• herding_gain: positive scalar controlling herding strength (default: 1.0)
• k_herding: alternative parameter name for herding_gain (for backward compatibility)

Returns

• ControlledTrajectorySheaf{T} representing the herding vehicle platoon system

Notes

The herding behavior emerges from the sheaf structure:

• When herding_gain = 0, vehicles evolve independently (decoupled)
• When herding_gain > 0, vehicles are attracted to each other with force proportional to their separation
• The sheaf Laplacian automatically generates the appropriate coupling terms

expected_feasible_dimension(ts::ControlledTrajectorySheaf{T};
                            rtol::Real=1e-10) where T
    -> Int

Return the expected dimension of the endpoint-preserving feasible trajectory family for ts, computed as

\[k m - \operatorname{rank}(\mathcal{C}_k),\]

where m is the control dimension and \mathcal{C}_k is the finite-horizon controllability matrix from finite_horizon_controllability.

feasible_control_trajectory_basis(ts::ControlledTrajectorySheaf{T},
                                  x1::AbstractVector{T},
                                  xk1::AbstractVector{T}) where T
    -> (z_p::Vector{T}, null_basis::Matrix{T})

Return one feasible trajectory hitting the endpoint states together with a basis matrix whose columns span all feasible perturbations that preserve those endpoint states.

The boundary conditions fix vertex 1 to the initial state x1 and vertex k+2 to the terminal state xk1 of the inner sheaf. The particular solution and each null-basis column are then converted to the public trajectory coordinates

\[z = (x_1, x_2, \dots, x_{k+1}, u_1, u_2, \dots, u_k).\]

The complete feasible set is { z_p + null_basis * α : α ∈ Rᵈ }.

Arguments

• ts — a ControlledTrajectorySheaf{T}.
• x1 — initial state vector of length ts.state_dim.
• xk1 — terminal state vector of length ts.state_dim.

Throws

• ArgumentError if length(x1) ≠ ts.state_dim.
• ArgumentError if length(xk1) ≠ ts.state_dim.

finite_horizon_controllability(ts::ControlledTrajectorySheaf{T}) where T
    -> Matrix{T}

Return the finite-horizon controllability matrix

\[\mathcal{C}_k = egin{bmatrix} A_d^{k-1} B_d & A_d^{k-2} B_d & \cdots & B_d nd{bmatrix}\]

for the k-step sampled dynamics stored in ts.

lqr_objective(ts::ControlledTrajectorySheaf{T},
              Q::AbstractMatrix{T},
              Ru::AbstractMatrix{T};
              Qf::Union{Nothing,AbstractMatrix{T}}=nothing,
              x_ref::Union{Nothing,AbstractMatrix{T}}=nothing,
              u_ref::Union{Nothing,AbstractMatrix{T}}=nothing) where T
    -> (H::SparseMatrixCSC{T,Int}, f::Vector{T}, c::T)

Assemble the quadratic objective

\[\tfrac{1}{2} z^\top H z + f^\top z + c\]

for the controlled sampled trajectory

\[z = (x_1, x_2, \dots, x_{k+1}, u_1, u_2, \dots, u_k).\]

The stagewise cost is

\[J(z) = \tfrac{1}{2} \sum_{t=1}^{k} \Bigl((x_t - \bar{x}_t)^\top Q (x_t - \bar{x}_t) + (u_t - \bar{u}_t)^\top R_u (u_t - \bar{u}_t)\Bigr) + \tfrac{1}{2} (x_{k+1} - \bar{x}_{k+1})^\top Q_f (x_{k+1} - \bar{x}_{k+1}).\]

The assembled H is sparse and block-diagonal: k copies of Q on the running state blocks, Qf on the terminal state block, and k copies of Ru on the control blocks.

Arguments

• ts — a ControlledTrajectorySheaf{T}.
• Q — $n \times n$ running state weight, positive semidefinite.
• Ru — $m \times m$ control weight, symmetric positive definite.
• Qf — $n \times n$ terminal state weight (defaults to Q when nothing).
• x_ref — $n \times (k+1)$ reference state trajectory (defaults to zeros).
• u_ref — $m \times k$ reference control trajectory (defaults to zeros).

Throws

• ArgumentError if Q is not $n \times n$.
• ArgumentError if Qf (when provided) is not $n \times n$.
• ArgumentError if Ru is not $m \times m$ or is not symmetric positive definite.
• ArgumentError if x_ref (when provided) is not $n \times (k+1)$.
• ArgumentError if u_ref (when provided) is not $m \times k$.

optimal_control_trajectory(ts::ControlledTrajectorySheaf{T},
                           x1::AbstractVector{T},
                           xk1::AbstractVector{T},
                           H::AbstractMatrix{T},
                           f::AbstractVector{T}=zeros(T, size(H, 1))) where T
    -> (z_opt::BlockVector, α_opt::Vector{T}, z_p::BlockVector, null_basis::AbstractMatrix{T})

Compute the minimum-cost feasible trajectory for the convex quadratic objective

\[\tfrac{1}{2} z^\top H z + f^\top z\]

subject to fixed initial state x1 and terminal state xk1.

The feasible set is the affine space $\{z_p + N\alpha : \alpha \in \mathbb{R}^r\}$ from feasible_control_trajectory_basis, and the optimizer satisfies the reduced first-order conditions

\[(N^\top H N)\, \alpha^\star = -N^\top (H z_p + f).\]

When the reduced Hessian $N^\top H N$ is singular, the minimum-norm optimizer is returned (using the pseudoinverse). If the reduced problem is unbounded below (i.e. the right-hand side has a component outside the range of the reduced Hessian), or if the reduced Hessian is not positive semidefinite, an ArgumentError is thrown.

Arguments

• ts — a ControlledTrajectorySheaf{T}.
• x1 — initial state vector of length ts.state_dim.
• xk1 — terminal state vector of length ts.state_dim.
• H — $p \times p$ positive-semidefinite Hessian matrix, where $p = (k+1) \cdot n + k \cdot m$.
• f — $p$-vector of linear objective coefficients (defaults to zero).

Returns

• z_opt — optimal trajectory as a BlockVector with k+1 state blocks followed by k control blocks.
• α_opt — optimal reduced coordinate vector.
• z_p — particular feasible trajectory as a BlockVector.
• null_basis — null-basis matrix N (columns span endpoint-preserving perturbations).

Throws

• ArgumentError if size(H) is not p × p.
• ArgumentError if length(f) is not p.
• ArgumentError if the reduced quadratic problem is unbounded below.
• ArgumentError if the reduced Hessian is not positive semidefinite.