Closed-loop MPC for multi-agent tracking sheaf problems by Optimax14 · Pull Request #65 · AlgebraicJulia/CellularSheaves.jl

Optimax14 · 2026-06-05T05:00:30Z

This PR implements a receding-horizon MPC algorithm for multi-agent tracking sheaf problems.

Summary

Added window_problem to take the entire horizon and slice it up into sub problems of size k.
Added pre-computation of sheaf Laplacian factorization, which turns each MPC step from resolving the open loop problem on a smaller horizon to a simple solve of the harmonic extension.
Added run_mpc_scenario with :cached (reuses the factorization + orthogonal projection) & :naive (recomputes the open loop problem for each MPC step)
Added ldlt_pinv_solve helpers to EuclideanSheaves and switched solve_quadratic_on_basis to use svd\

Tests

40 new tests in test/network_sheaves/MultiAgentTrackingMPC.jl.

Docs

New literate page in docs/literature/control/mpc_target_tracking.jl showcasing the results of this PR and the application of the algorithm in different scenarios.

Closes #62

…ion, run_mpc_scenario

codecov · 2026-06-05T05:09:19Z

Codecov Report

❌ Patch coverage is 98.87006% with 2 lines in your changes missing coverage. Please review.
✅ Project coverage is 75.06%. Comparing base (c379f01) to head (a2ffcf4).

Files with missing lines	Patch %	Lines
src/ControlSheaves/MultiAgentTracking.jl	99.24%	1 Missing ⚠️
src/network_sheaves/EuclideanSheaves.jl	97.50%	1 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main      #65      +/-   ##
==========================================
+ Coverage   73.60%   75.06%   +1.45%     
==========================================
  Files          36       36              
  Lines        2550     2707     +157     
==========================================
+ Hits         1877     2032     +155     
- Misses        673      675       +2

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

Copilot

Copilot encountered an error and was unable to review this pull request. You can try again by re-requesting a review.

Copilot

Pull request overview

Copilot reviewed 8 out of 8 changed files in this pull request and generated 5 comments.

+    for i in 1:prob.n_agents
+        va = agent_vertex(prob, i, 0)
+        for c in 1:size(prob.Ad[i], 1)
+            pinned[offsets[va] + c] = x_agents[i][c]
+        end
+    end


+    for j in 1:prob.n_targets, ts in 0:prob.k
+        vt  = target_vertex(prob, j, ts)
+        val = target_window[j][ts + 1]
+        for c in 1:length(val)
+            pinned[offsets[vt] + c] = val[c]
+        end
+    end


+    nx          = [size(prob.Ad[i], 1) for i in 1:prob.n_agents]
+    x_now       = copy.(x0_agents)
+    agent_trajs = [Matrix{Float64}(undef, k+1, nx[i]) for i in 1:prob.n_agents]
+    for i in 1:prob.n_agents; agent_trajs[i][1, :] = x_now[i]; end


+            x_B = Vector{Float64}(undef, length(op.boundary))
+            p = 1
+            for i in 1:prob.n_agents, c in 1:nx[i]; x_B[p] = x_now[i][c]; p += 1; end
+            for ts in 0:op.window, j in 1:prob.n_targets
+                for c in eachindex(local_targets[j][ts+1]); x_B[p] = local_targets[j][ts+1][c]; p += 1; end
+            end
+            z        = op(x_B)


jpfairbanks · 2026-06-05T18:11:29Z

-    α_opt = pinv(QN) * rhs
+    QN = Matrix(N' * Q * N)
+    rhs = -(N' * Q * point)
+    α_opt = svd(QN) \ rhs


Why is this using svd instead of ldl_nullspace?

jpfairbanks · 2026-06-05T18:42:37Z

+# Overload that pins individual stalk components rather than whole vertices.
+# Works directly on the flat Laplacian since the vertex-level helpers above
+# cannot express partial pinning.
+function harmonic_extension(s::EuclideanSheaf, boundary::Dict{Int,Float64})


Is this dead code? I thought we discussed a way to avoid the need for partial pinning that breaks the abstraction of harmonic extension.

jpfairbanks · 2026-06-05T18:45:51Z

+        _default_ldlt_nullity_threshold(maximum(i -> abs(D[i, i]), 1:n; init=0.0)) : tol
+    c = F.P' \ b
+    z = F.L \ c
+    w = ldiv_diag_pinv(D, z, threshold)


this should use the same ldiv_diag_pinv! as above. If we are reusing the F factorization object, then we might as well bundle a w vector buffer to reuse as well.

jpfairbanks · 2026-06-05T18:47:26Z

+    n = size(D, 1)
+    threshold = isnothing(tol) ?
+        _default_ldlt_nullity_threshold(maximum(i -> abs(D[i, i]), 1:n; init=0.0)) : tol
+    c = F.P' \ b


you can do an in place ldiv and reuse all the memory. 0 allocations per step would be ideal

jpfairbanks · 2026-06-05T18:47:42Z

+    c = F.P' \ b
+    z = F.L \ c
+    w = ldiv_diag_pinv(D, z, threshold)
+    return F.P \ (F.L' \ w)


in place ldiv!

jpfairbanks · 2026-06-05T18:56:11Z

+    L_II = sparse(L[interior, interior])
+    L_IB = sparse(L[interior, boundary])
+    F    = ldlt!(ChordalLDLt(L_II), RowMaximum(); check=false)
+    tol  = sqrt(eps(Float64)) * max(1.0, maximum(i -> abs(F.D[i,i]), 1:size(F.D,1); init=0.0))


I think this logic is repeated at least one other place and so should be a function.

jpfairbanks · 2026-06-05T18:57:07Z

+        Fq     = ldlt!(DenseLDLtPivoted(NtQ_II * N), RowMaximum(); check=false)
+        rhs    = NtQ_II * H
+        coeff  = similar(rhs)
+        for col in axes(rhs, 2); coeff[:, col] = ldlt_pinv_solve(Fq, rhs[:, col]); end


Is this computing a nullspace with multiple right hand sides? Maybe this should be its own method?

jpfairbanks · 2026-06-05T19:00:42Z

+    return TrackingExtension(M, boundary, interior, stalks, sparse(L), size(N, 2), W)
+end
+
+function (op::TrackingExtension)(x_B::AbstractVector)


nice use of function call overloading.

jpfairbanks · 2026-06-05T19:01:00Z

+
+function (op::TrackingExtension)(x_B::AbstractVector)
+    @argcheck length(x_B) == length(op.boundary)
+    z = zeros(sum(op.stalks))


This could be cached in the cache and returned by reference.

or use a mul!(z::AbstractVector, op::TrackingExtension, x_B::AbstractVector)

jpfairbanks · 2026-06-05T19:02:42Z

+            L        = sheaf_laplacian_matrix_direct(sheaf)
+            t == 0 && (nd = size(N, 2))
+            residual = sqrt(max(0.0, dot(z_arr, L * z_arr)))
+            x_now    = _apply_first_control(local_prob, z_arr, sheaf.vertex_stalks)


try to minimize the repeated logic across the use_cached switch

Optimax14 added 4 commits June 5, 2026 00:04

add ldlt_pinv_solve helpers and use svd\ in solve_quadratic_on_basis

0d6d41a

add MPC support to MultiAgentTracking: window_problem, TrackingExtens…

8de7ee0

…ion, run_mpc_scenario

add MPC tests for MultiAgentTracking

f14ba58

add MPC literate doc and register in docs/make.jl

8990b4d

Optimax14 force-pushed the itay/mpc-multi-agent branch from 61d075f to 8990b4d Compare June 5, 2026 05:05

fix p_z plot not rendering in docs

a2ffcf4

jpfairbanks requested a review from Copilot June 5, 2026 18:15

Copilot started reviewing on behalf of jpfairbanks June 5, 2026 18:16 View session

Copilot AI reviewed Jun 5, 2026

View reviewed changes

jpfairbanks requested a review from Copilot June 5, 2026 18:40

Copilot started reviewing on behalf of jpfairbanks June 5, 2026 18:41 View session

Copilot AI reviewed Jun 5, 2026

View reviewed changes

jpfairbanks requested changes Jun 5, 2026

View reviewed changes

Optimax14 self-assigned this Jun 5, 2026

Conversation

Optimax14 commented Jun 5, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Summary

Tests

Docs

Uh oh!

codecov Bot commented Jun 5, 2026

Codecov Report

Uh oh!

Copilot AI left a comment

Choose a reason for hiding this comment

Uh oh!

Copilot AI left a comment

Choose a reason for hiding this comment

Pull request overview

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

Optimax14 commented Jun 5, 2026 •

edited

Loading