new routing algorithm

@miSim/initialize.m

@@ -1,4 +1,4 @@
-function [obj] = initialize(obj, domain, agents, barrierGain, barrierExponent, minAlt, timestep, maxIter, obstacles, makePlots, makeVideo, useDoubleIntegrator, dampingCoeff, useFixedTopology, optimizeSensorPointing, useSinrComms, sinrThreshold, pathLossExponent, ambientTemp, centerFreq, bandwidth, useRoutingTopology, routingFlowThreshold)
+function [obj] = initialize(obj, domain, agents, barrierGain, barrierExponent, minAlt, timestep, maxIter, obstacles, makePlots, makeVideo, useDoubleIntegrator, dampingCoeff, useFixedTopology, optimizeSensorPointing, useSinrComms, sinrThreshold, pathLossExponent, ambientTemp, centerFreq, bandwidth, useRoutingTopology)
     arguments (Input)
         obj (1, 1) {mustBeA(obj, "miSim")};
         domain (1, 1) {mustBeGeometry};
@@ -22,7 +22,6 @@ function [obj] = initialize(obj, domain, agents, barrierGain, barrierExponent, m
         centerFreq (1, 1) double = 2.4e9;
         bandwidth (1, 1) double = 20e6;
         useRoutingTopology (1, 1) logical = false;
-        routingFlowThreshold (1, 1) double = 0.02;
     end
     arguments (Output)
         obj (1, 1) {mustBeA(obj, "miSim")};
@@ -114,23 +113,49 @@ function [obj] = initialize(obj, domain, agents, barrierGain, barrierExponent, m
     obj.centerFreq = centerFreq;
     obj.bandwidth = bandwidth;
 
-    % Set topology-selection algorithm and routing parameters
+    % Set topology-selection algorithm; fail fast on invalid mode combinations
+    assert(~(useRoutingTopology && useFixedTopology), ...
+        "useRoutingTopology and useFixedTopology are mutually exclusive topology selectors");
+    assert(~useRoutingTopology || useSinrComms, ...
+        "useRoutingTopology requires the SINR comms model (useSinrComms)");
     obj.useRoutingTopology = useRoutingTopology;
-    obj.routingFlowThreshold = routingFlowThreshold;
 
-    % Compute adjacency matrix and network topology. Lesser-neighbor always
-    % provides the low-level connectivity topology (basic traffic must be
-    % routable between any two drones at all times); routing mode overlays
-    % the bulk-data links on top of it. The routing LP uses linprog (no
-    % codegen support); the compiled path never sets the flag, so linprog
-    % is never compiled.
+    % Compute adjacency matrix and network topology. The Capacity-Aware
+    % Lesser Sink Neighbor algorithm replaces lesser-neighbor when selected
+    % (MATLAB sim only; the compiled path never sets the flag).
     obj = obj.updateAdjacency();
     if obj.useFixedTopology
         obj.constraintAdjacencyMatrix = obj.adjacency;
+    elseif coder.target('MATLAB') && obj.useRoutingTopology
+        obj = obj.lesserSinkNeighbor();
     else
         obj = obj.lesserNeighbor();
-        if coder.target('MATLAB') && obj.useRoutingTopology
-            obj = obj.routeTopology();
+    end
+
+    % Guard against the SINR CBF deadlock band. Connectivity admits links
+    % with SINR anywhere above sinrThreshold, but constrainMotion maintains
+    % every selected link at sinrThreshold + 3 dB (paddingDb — keep in
+    % sync) in BOTH directions. A link selected inside that 3 dB band is
+    % born with a violated barrier; the QP can then go infeasible and hold
+    % positions, and since a frozen geometry re-selects the same link, the
+    % fleet deadlocks permanently. This catches the initial configuration;
+    % geometry evolution can still enter the band mid-run (fixing that
+    % requires pad-aware selection or a bounded-recovery barrier floor).
+    if coder.target('MATLAB') && obj.useSinrComms
+        gammaCbfLin = 10^((obj.sinrThreshold + 3.0) / 10);
+        posChk = zeros(size(obj.agents, 1), 3);
+        for kk = 1:size(obj.agents, 1)
+            posChk(kk, :) = obj.agents{kk}.pos;
+        end
+        for ii = 1:(size(obj.agents, 1) - 1)
+            for jj = (ii + 1):size(obj.agents, 1)
+                if obj.constraintAdjacencyMatrix(ii, jj)
+                    sinrBinding = min(obj.sinrLink(posChk, ii, jj), obj.sinrLink(posChk, jj, ii));
+                    assert(sinrBinding >= gammaCbfLin, ...
+                        "Initial topology link %d-%d has SINR %.1f dB, inside the CBF maintenance band (sinrThreshold %.1f dB + 3 dB padding): the QP starts with a violated comms barrier and can deadlock. Lower sinrThreshold or adjust the initial geometry.", ...
+                        ii, jj, 10 * log10(sinrBinding), obj.sinrThreshold);
+                end
+            end
         end
     end
 

@miSim/initializeFromCsv.m

@@ -147,17 +147,12 @@ else
     BANDWIDTH = 20e6;
 end
 
-% Routing topology selection (optional columns; MATLAB sim only)
+% Routing topology selection (optional column; MATLAB sim only)
 if isfield(scenario, 'useRoutingTopology')
     USE_ROUTING_TOPOLOGY = logical(scenario.useRoutingTopology);
 else
     USE_ROUTING_TOPOLOGY = false;
 end
-if isfield(scenario, 'routingFlowThreshold')
-    ROUTING_FLOW_THRESHOLD = scenario.routingFlowThreshold;
-else
-    ROUTING_FLOW_THRESHOLD = 0.02;
-end
 
 % ---- Build domain --------------------------------------------------------
 dom = rectangularPrism;
@@ -208,6 +203,6 @@ obj = obj.initialize(dom, agentList, BARRIER_GAIN, BARRIER_EXPONENT, ...
                      MIN_ALT, TIMESTEP, MAX_ITER, obstacleList, false, false, ...
                      USE_DOUBLE_INTEGRATOR, DAMPING_COEFF, USE_FIXED_TOPOLOGY, false, ...
                      USE_SINR_COMMS, SINR_THRESHOLD, PATH_LOSS_EXPONENT, AMBIENT_TEMP, CENTER_FREQ, BANDWIDTH, ...
-                     USE_ROUTING_TOPOLOGY, ROUTING_FLOW_THRESHOLD);
+                     USE_ROUTING_TOPOLOGY);
 
 end

@miSim/initializeFromInits.m

@@ -116,17 +116,12 @@ if isfield(inits, 'useRoutingTopology')
 else
     useRoutingTopology = false;
 end
-if isfield(inits, 'routingFlowThreshold')
-    routingFlowThreshold = inits.routingFlowThreshold;
-else
-    routingFlowThreshold = 0.02;
-end
 
 % ---- Initialize simulation (plots and video disabled) ------------------------
 obj = obj.initialize(dom, agentList, inits.barrierGain, inits.barrierExponent, ...
                      inits.minAlt, inits.timestep, inits.maxIter, obstacleList, ...
                      false, false, useDoubleIntegrator, dampingCoeff, useFixedTopology, false, ...
                      useSinrComms, sinrThreshold, pathLossExponent, ambientTemp, centerFreq, bandwidth, ...
-                     useRoutingTopology, routingFlowThreshold);
+                     useRoutingTopology);
 
 end

@miSim/lesserSinkNeighbor.m

@@ -0,0 +1,372 @@
+function obj = lesserSinkNeighbor(obj)
+% LESSERSINKNEIGHBOR  Capacity-Aware Lesser Sink Neighbor algorithm (CALSN).
+%
+% Implements the algorithm specified in new/new_routing_algorithm.pdf as the
+% selectable alternative to lesserNeighbor. UAVs 1 and N are base stations;
+% every other UAV generates r = 1 unit of data per step, and all data must
+% reach a base station. The algorithm outputs both a sparse maintained
+% topology and a routing table.
+%
+% 1. Base-station potentials phi via distributed Bellman-Ford over the
+%    current communication neighborhoods (obj.adjacency), with link cost
+%    c_ij = 1/(C_ij + eps), where C_ij = bandwidth * log2(1 + SINR_ij) is
+%    the Shannon capacity of direction i -> j.
+% 2. Lesser sink neighbors L_i replace LNA's "j < i" rule:
+%        L_i = { j in N_i : phi_j < phi_i, or phi_j == phi_i and j < i }.
+%    The equal-potential lexicographic tie-break is a documented deviation
+%    from the spec's strict inequality: at an exact potential watershed
+%    (e.g. a symmetric chain) the strict rule would split A_control into
+%    one tree per base station; the tie-break keeps the airborne graph
+%    connected, which the mission requires.
+% 3. LNA-style sparsity: one parent per connected component of the
+%    L_i-induced subgraph (only links with positive capacity in at least
+%    one direction count as edges), chosen by
+%        p = argmin_{j in C} ( phi_j + 1/(C_ij + eps) ).
+%    A_control(i,p) = A_control(p,i) = 1 (diagonal = I).
+%    BASIN EXTENSION (documented deviation): with two potential minima
+%    (the two bases), one-parent-per-component alone yields a spanning
+%    forest with up to two disjoint base-rooted trees whenever a UAV's
+%    lesser sink neighbors collapse into a single component (dense usable
+%    graphs). To restore LNA's global-connectivity property, each UAV
+%    selects one representative per component AND per base-station BASIN
+%    (the base a neighbor's potential descent terminates at). Interior
+%    UAVs are unaffected (their whole neighborhood shares one basin);
+%    watershed UAVs whose neighborhoods span both basins gain a second
+%    parent, which is exactly the bridge between the two trees. A final
+%    stitch guard covers degenerate geometries: if A_control still has
+%    multiple components, the lowest-cost feasible inter-component link
+%    is added (maintained for connectivity only; carries no flow).
+% 4. Routing, children first (decreasing phi; equal phi processed in
+%    decreasing index order so tie-break children precede their parents):
+%    aggregate demand D_i = r + inflows, split capacity-proportionally
+%    over parents; saturated links are capped at capacity and the residual
+%    redistributed over unsaturated parents; any demand beyond the total
+%    parent capacity is reported as unrouted and the UAV marked infeasible.
+%
+% Outputs (set on obj):
+%   constraintAdjacencyMatrix  A_control: undirected maintained topology
+%   routingFlows               F: flow on each directed link (units of r)
+%   routingFractions           R = F/D: share of the sender's demand, 0..1
+%   routingAdjacencyMatrix     F > 0 (directed, child -> parent)
+%   routingPotentials          phi (base-station potential per UAV)
+%   routingUnrouted            unrouted demand per UAV (zeros when feasible)
+
+    arguments (Input)
+        obj (1, 1) {mustBeA(obj, "miSim")};
+    end
+    arguments (Output)
+        obj (1, 1) {mustBeA(obj, "miSim")};
+    end
+
+    nAgents = size(obj.agents, 1);
+    r = 1.0;         % data generated per non-base UAV per step (abstract unit)
+    EPS_CAP = 1e-12; % epsilon in c = 1/(C + eps); avoids division by zero
+    TOL = 1e-12;     % residual-demand tolerance for the saturation loop
+
+    isBase = false(nAgents, 1);
+    isBase(1) = true;
+    isBase(nAgents) = true; % base stations are UAVs 1 and N
+
+    positions = zeros(nAgents, 3);
+    for kk = 1:nAgents
+        positions(kk, :) = obj.agents{kk}.pos;
+    end
+
+    % Directed Shannon link capacities over the current communication
+    % neighborhoods: C(i, j) = capacity of sending from i to j (receiver j)
+    C = zeros(nAgents);
+    for ii = 1:nAgents
+        for jj = 1:nAgents
+            if jj ~= ii && obj.adjacency(ii, jj)
+                C(ii, jj) = obj.bandwidth * log2(1 + obj.sinrLink(positions, jj, ii));
+            end
+        end
+    end
+
+    % ---- Base-station potentials (Bellman-Ford, at most N-1 sweeps) -----
+    phi = inf(nAgents, 1);
+    phi(isBase) = 0;
+    for sweep = 1:max(nAgents - 1, 1)
+        changed = false;
+        for ii = 1:nAgents
+            if isBase(ii)
+                continue;
+            end
+            newPhi = inf;
+            for jj = 1:nAgents
+                if jj == ii || ~obj.adjacency(ii, jj)
+                    continue;
+                end
+                newPhi = min(newPhi, phi(jj) + 1 / (C(ii, jj) + EPS_CAP));
+            end
+            if newPhi ~= phi(ii)
+                phi(ii) = newPhi;
+                changed = true;
+            end
+        end
+        if ~changed
+            break;
+        end
+    end
+
+    % ---- Base-station basins (for the basin extension) -------------------
+    % Each finite-potential UAV descends its argmin neighbor chain; phi is
+    % strictly decreasing along it, so the chain terminates at a base. The
+    % basin label records which base that is.
+    descent = zeros(nAgents, 1);
+    for ii = 1:nAgents
+        if isBase(ii) || isinf(phi(ii))
+            continue;
+        end
+        bestMetric = inf;
+        for jj = 1:nAgents
+            if jj == ii || ~obj.adjacency(ii, jj)
+                continue;
+            end
+            metric = phi(jj) + 1 / (C(ii, jj) + EPS_CAP);
+            if metric < bestMetric
+                bestMetric = metric;
+                descent(ii) = jj;
+            end
+        end
+    end
+    basin = zeros(nAgents, 1);
+    basin(1) = 1;
+    basin(nAgents) = nAgents;
+    for ii = 1:nAgents
+        if isBase(ii) || isinf(phi(ii))
+            continue;
+        end
+        kk = ii;
+        while ~isBase(kk)
+            kk = descent(kk);
+        end
+        basin(ii) = kk;
+    end
+
+    % ---- Parent selection (one per lesser-sink component and basin) ------
+    Acontrol = logical(eye(nAgents));
+    isParent = false(nAgents); % isParent(i, p): p is a selected parent of i
+    for ii = 1:nAgents
+        if isBase(ii) || isinf(phi(ii))
+            % Base stations select no parents; an infinite potential means
+            % no current local route toward a base station
+            continue;
+        end
+
+        % Lesser sink neighbors (with the equal-potential tie-break)
+        lsnBuf = zeros(1, nAgents);
+        lsnCount = 0;
+        for jj = 1:nAgents
+            if jj == ii || ~obj.adjacency(ii, jj)
+                continue;
+            end
+            if phi(jj) < phi(ii) || (phi(jj) == phi(ii) && jj < ii)
+                lsnCount = lsnCount + 1;
+                lsnBuf(lsnCount) = jj;
+            end
+        end
+        if lsnCount == 0
+            continue;
+        end
+        lsn = lsnBuf(1:lsnCount);
+
+        % Subgraph induced by the lesser sink neighbors in the usable
+        % communication graph: only links with positive capacity in at
+        % least one direction count as edges (a zero-capacity link cannot
+        % carry data, so it must not merge two components)
+        subgraphAdjacency = false(lsnCount);
+        for aa = 1:lsnCount
+            for bb = (aa + 1):lsnCount
+                ja = lsn(aa);
+                jb = lsn(bb);
+                if obj.adjacency(ja, jb) && (C(ja, jb) > 0 || C(jb, ja) > 0)
+                    subgraphAdjacency(aa, bb) = true;
+                    subgraphAdjacency(bb, aa) = true;
+                end
+            end
+        end
+
+        % One parent per connected component and per base-station basin
+        % present in it (basin extension), selected by capacity and
+        % potential: p = argmin( phi_j + 1/(C_ij + eps) )
+        visited = false(1, lsnCount);
+        for aa = 1:lsnCount
+            if visited(aa)
+                continue;
+            end
+            component = bfs(subgraphAdjacency, aa);
+            visited(component) = true;
+
+            for bb = [1, nAgents] % the two possible basin labels
+                bestMetric = inf;
+                bestParent = 0;
+                for cc = 1:numel(component)
+                    jj = lsn(component(cc));
+                    if basin(jj) ~= bb
+                        continue;
+                    end
+                    metric = phi(jj) + 1 / (C(ii, jj) + EPS_CAP);
+                    if metric < bestMetric
+                        bestMetric = metric;
+                        bestParent = jj;
+                    end
+                end
+                if bestParent > 0
+                    isParent(ii, bestParent) = true;
+                    Acontrol(ii, bestParent) = true;
+                    Acontrol(bestParent, ii) = true;
+                end
+            end
+        end
+    end
+
+    % ---- Connectivity stitch guard ---------------------------------------
+    % The basin extension bridges the two base-rooted trees through
+    % watershed UAVs whenever each base borders a routing basin. Degenerate
+    % geometries (e.g. a base whose only feasible neighbor is the other
+    % base) can still leave A_control split; add the lowest-cost feasible
+    % link between components until one remains. Stitched links are
+    % maintained for connectivity only and carry no routed flow. If the
+    % feasibility graph itself is disconnected no stitch exists; validate()
+    % reports that case.
+    while true
+        comp = zeros(1, nAgents);
+        nComp = 0;
+        for ii = 1:nAgents
+            if comp(ii) == 0
+                nComp = nComp + 1;
+                comp(bfs(Acontrol, ii)) = nComp;
+            end
+        end
+        if nComp <= 1
+            break;
+        end
+        bestMetric = inf;
+        bestU = 0;
+        bestV = 0;
+        for uu = 1:nAgents
+            for vv = (uu + 1):nAgents
+                if comp(uu) ~= comp(vv) && obj.adjacency(uu, vv)
+                    % Bidirectional link cost (finite even for infinite
+                    % potentials, so isolated islands remain stitchable)
+                    metric = 1 / (C(uu, vv) + EPS_CAP) + 1 / (C(vv, uu) + EPS_CAP);
+                    if metric < bestMetric
+                        bestMetric = metric;
+                        bestU = uu;
+                        bestV = vv;
+                    end
+                end
+            end
+        end
+        if bestU == 0
+            break; % feasibility graph disconnected: nothing to stitch with
+        end
+        Acontrol(bestU, bestV) = true;
+        Acontrol(bestV, bestU) = true;
+    end
+
+    % ---- Routing (children first: decreasing potential order) -----------
+    % Equal potentials are processed in decreasing index order: with the
+    % lexicographic tie-break, flow within an equal-potential group goes
+    % from higher to lower index, so higher indices must aggregate first.
+    F = zeros(nAgents);
+    R = zeros(nAgents);
+    unrouted = zeros(nAgents, 1);
+
+    nonBase = find(~isBase);
+    [~, ord] = sortrows([phi(nonBase), nonBase], [-1, -2]);
+    order = nonBase(ord);
+
+    for oo = 1:numel(order)
+        ii = order(oo);
+
+        % Aggregated outgoing demand: own data plus everything received
+        D = r + sum(F(:, ii));
+        parents = find(isParent(ii, :));
+        if D == 0 || isempty(parents)
+            continue;
+        end
+
+        % Capacity-proportional split, saturating full links first and
+        % redistributing the remainder across parents with remaining
+        % capacity. Each pass either routes all residual demand or
+        % saturates at least one parent, so the loop terminates.
+        capP = C(ii, parents);
+        flows = zeros(1, numel(parents));
+        active = true(1, numel(parents));
+        residual = D;
+        while residual > TOL && any(active)
+            totalCap = sum(capP(active));
+            if totalCap <= 0
+                break;
+            end
+            alloc = zeros(1, numel(parents));
+            alloc(active) = residual .* capP(active) ./ totalCap;
+            flows = flows + alloc;
+            saturated = active & (flows > capP);
+            if any(saturated)
+                residual = sum(flows(saturated) - capP(saturated));
+                flows(saturated) = capP(saturated);
+                active(saturated) = false;
+            else
+                residual = 0;
+            end
+        end
+        if residual > TOL
+            % Demand exceeds the total parent capacity: route as much as
+            % possible and report the remainder as unrouted (infeasible)
+            unrouted(ii) = residual;
+            warning("miSim:routingInfeasible", ...
+                "UAV %d demand %.3g exceeds total parent capacity; %.3g units unrouted", ...
+                ii, D, residual);
+        end
+
+        F(ii, parents) = flows;
+        R(ii, parents) = flows ./ D;
+    end
+
+    % ---- Outputs ---------------------------------------------------------
+    obj.constraintAdjacencyMatrix = Acontrol;
+    obj.routingFlows = F;
+    obj.routingFractions = R;
+    obj.routingAdjacencyMatrix = F > 0;
+    obj.routingPotentials = phi;
+    obj.routingUnrouted = unrouted;
+end
+
+function cComp = bfs(subgraphAdjacency, startIdx)
+    % Breadth-first search returning the connected component of startIdx.
+    % Pre-allocated queue and component buffer with head/tail pointers to
+    % avoid element deletion and dynamic array growth (same pattern as the
+    % local bfs in lesserNeighbor.m).
+    n = size(subgraphAdjacency, 1);
+    visited = false(1, n);
+
+    queue    = zeros(1, n);
+    cCompBuf = zeros(1, n);
+    qHead = 1;
+    qTail = 2;
+    queue(1)    = startIdx;
+    cCompBuf(1) = startIdx;
+    cSize = 1;
+    visited(startIdx) = true;
+
+    while qHead < qTail
+        current = queue(qHead);
+        qHead = qHead + 1;
+
+        neighbors = find(subgraphAdjacency(current, :));
+        for kk = 1:numel(neighbors)
+            neighbor = neighbors(kk);
+            if ~visited(neighbor)
+                visited(neighbor) = true;
+                cCompBuf(cSize + 1) = neighbor;
+                cSize = cSize + 1;
+                queue(qTail) = neighbor;
+                qTail = qTail + 1;
+            end
+        end
+    end
+    cComp = sort(cCompBuf(1:cSize));
+end

@miSim/miSim.m

@@ -31,13 +31,17 @@ classdef miSim
         centerFreq = 2.4e9;       % carrier frequency (Hz); sets path-loss offset K = (4*pi*f_c/c)^2
         bandwidth = 20e6;         % channel bandwidth (Hz) for thermal noise N = k_B*T*B
         % Topology-selection algorithm (requires useSinrComms; MATLAB sim only).
-        % false -> lesser-neighbor; true -> min-max-workload LP routing, which
-        % picks the links that route all sensing data to the network endpoints
-        % (agents 1 and the largest odd index) with minimal peak node workload.
+        % false -> lesser-neighbor; true -> Capacity-Aware Lesser Sink Neighbor
+        % (CALSN, see new/new_routing_algorithm.pdf): base-station potentials
+        % replace the index-based "lesser" rule, and the maintained topology
+        % comes with a routing table (flows F, fractions R) that sends every
+        % UAV's data toward base stations 1 and N.
         useRoutingTopology = false;
-        routingFlowThreshold = 0.02;  % prune links carrying less flow than this (units of one node's production)
-        routingAdjacencyMatrix = false(0, 0); % directed adjacency chosen by the routing LP: (i, j) = link i -> j
-        routingFlows = zeros(0, 0);   % data flow carried by each directed link (units/timestep)
+        routingAdjacencyMatrix = false(0, 0); % directed flow links (i, j) = child i sends to parent j (F > 0)
+        routingFlows = zeros(0, 0);      % F: flow on each directed link (units of r = 1 per UAV per step)
+        routingFractions = zeros(0, 0);  % R = F/D: share of the sender's aggregated demand per link (0..1)
+        routingPotentials = zeros(0, 1); % phi: Bellman-Ford base-station potential per UAV
+        routingUnrouted = zeros(0, 1);   % unrouted demand per UAV (zeros when routing is feasible)
         artifactName = "";
         f; % main plotting tiled layout figure
         fPerf; % performance plot figure
@@ -87,13 +91,13 @@ classdef miSim
             obj.obstacles = {rectangularPrism};
             obj.agents = {agent};
         end
-        [obj] = initialize(obj, domain, agents, barrierGain, barrierExponent, minAlt, timestep, maxIter, obstacles, makePlots, makeVideo, useDoubleIntegrator, dampingCoeff, useFixedTopology, optimizeSensorPointing, useSinrComms, sinrThreshold, pathLossExponent, ambientTemp, centerFreq, bandwidth, useRoutingTopology, routingFlowThreshold);
+        [obj] = initialize(obj, domain, agents, barrierGain, barrierExponent, minAlt, timestep, maxIter, obstacles, makePlots, makeVideo, useDoubleIntegrator, dampingCoeff, useFixedTopology, optimizeSensorPointing, useSinrComms, sinrThreshold, pathLossExponent, ambientTemp, centerFreq, bandwidth, useRoutingTopology);
         [obj] = initializeFromCsv(obj, csvPath);
         [obj] = initializeFromInits(obj, initsPath);
         [obj] = plotFromSimHist(obj, initsPath, histPath);
         [obj] = run(obj);
         [obj] = lesserNeighbor(obj);
-        [obj] = routeTopology(obj);
+        [obj] = lesserSinkNeighbor(obj);
         [obj] = constrainMotion(obj);
         [sinr, grad] = sinrLink(obj, positions, rxIdx, txIdx);
         [obj] = partition(obj);

@miSim/plotGraph.m

@@ -6,11 +6,10 @@ function obj = plotGraph(obj)
         obj (1, 1) {mustBeA(obj, "miSim")};
     end
 
-    % The maintained (low-level connectivity) links are always drawn as a
-    % dotted undirected graph: basic traffic must be routable between any
-    % two drones at all times. In routing mode the bulk-data links selected
-    % by the LP are overlaid as SOLID DIRECTED edges (arrow = transmitter ->
-    % receiver) on the same node layout.
+    % The maintained links are always drawn as a dotted undirected graph.
+    % In routing (CALSN) mode the flow-carrying links are overlaid as SOLID
+    % DIRECTED edges (arrow = child -> parent, toward a base station) on
+    % the same node layout.
     G = graph(obj.constraintAdjacencyMatrix, "omitselfloops");
 
     % Resolve the target axes for the network graph tile
@@ -24,7 +23,7 @@ function obj = plotGraph(obj)
     hold(ax, "on");
     o = plot(ax, G, "LineStyle", "--", "EdgeColor", "g", "NodeColor", "k", "LineWidth", 2);
     if obj.useRoutingTopology
-        % Bulk overlay: reuse the base layout's node coordinates and hide
+        % Flow overlay: reuse the base layout's node coordinates and hide
         % the overlay's own nodes/labels so only its edges show. Sizeable
         % arrowheads so the flow direction reads at tile size.
         D = digraph(obj.routingAdjacencyMatrix);
@@ -51,11 +50,11 @@ function obj = plotGraph(obj)
             positions(kk, :) = obj.agents{kk}.pos;
         end
         % Labels are bare numbers (no units) so they stay compact; GraphPlot
-        % renders them horizontally, left-to-right. Dotted (low-level) edges
+        % renders them horizontally, left-to-right. Dotted (maintained) edges
         % are labelled with the binding (worse) directional SINR — a link is
         % feasible only if BOTH directions clear the threshold. Links that
-        % also carry a bulk overlay leave the label to the solid edge, which
-        % shows the SINR of the transmission direction.
+        % also carry flow leave the label to the solid directed edge, which
+        % shows the flow value (units of r) sent child -> parent.
         edgeLabels = strings(size(endNodes, 1), 1);
         for ee = 1:size(endNodes, 1)
             lo = endNodes(ee, 1);
@@ -70,13 +69,12 @@ function obj = plotGraph(obj)
         end
         o(1).EdgeLabel = edgeLabels;
         if obj.useRoutingTopology
-            endNodesD = D.Edges.EndNodes;   % each row [tx rx]
-            bulkLabels = strings(size(endNodesD, 1), 1);
+            endNodesD = D.Edges.EndNodes;   % each row [child parent]
+            flowLabels = strings(size(endNodesD, 1), 1);
             for ee = 1:size(endNodesD, 1)
-                sinrTx = obj.sinrLink(positions, endNodesD(ee, 2), endNodesD(ee, 1)); % rx, tx
-                bulkLabels(ee) = sprintf("%.1f", 10 * log10(sinrTx));
+                flowLabels(ee) = sprintf("%.2f", obj.routingFlows(endNodesD(ee, 1), endNodesD(ee, 2)));
             end
-            o(2).EdgeLabel = bulkLabels;
+            o(2).EdgeLabel = flowLabels;
         end
         % Threshold goes in the existing single-line title (not a subtitle) so
         % the small tile is not compressed further, and the hover toolbar is

@miSim/routeTopology.m

@@ -1,146 +0,0 @@
-function obj = routeTopology(obj)
-% ROUTETOPOLOGY  Select the network topology by solving a min-max-workload
-% routing LP (alternative to the lesser-neighbor algorithm).
-%
-% Every agent produces one unit of sensing data per timestep, all of which
-% must reach one of the network endpoints (agents 1 and the largest odd
-% index), which output data off-network for free (bottomless sinks). Data
-% may split across paths, but relayed data accumulates at intermediate
-% nodes, so each node experiences a workload: the share of a timestep it
-% spends transmitting and receiving, with per-link effort flow/capacity and
-% Shannon capacity c = log2(1 + SINR). The LP minimizes the maximum
-% workload over all nodes, since the busiest node bounds system capacity.
-%
-% Candidate links are the threshold-feasible pairs in obj.adjacency (both
-% directional SINRs clear sinrThreshold, per updateAdjacency). After the
-% first solve, links carrying negligible flow (< routingFlowThreshold) are
-% disallowed and the LP is re-solved so they cannot reappear: feasible but
-% impractical links add system complexity for very little gain.
-%
-% The routed links are a HIGH-QUALITY OVERLAY on top of the low-level
-% connectivity topology, not a replacement for it: basic command/telemetry
-% traffic must always be routable between any two drones, so the caller
-% runs lesserNeighbor first and this method ORs the bulk links into the
-% maintained set. Both tiers are held above threshold by the CBF.
-%
-% Outputs (set on obj):
-%   routingAdjacencyMatrix    directed adjacency, (i, j) = bulk link i -> j
-%   routingFlows              flow carried by each directed link
-%   constraintAdjacencyMatrix existing (lesser-neighbor) links OR the
-%                             symmetrized routed links (the CBF maintains
-%                             both directions of every link)
-%
-% Uses linprog (Optimization Toolbox): MATLAB simulation only, no codegen.
-
-    arguments (Input)
-        obj (1, 1) {mustBeA(obj, "miSim")};
-    end
-    arguments (Output)
-        obj (1, 1) {mustBeA(obj, "miSim")};
-    end
-
-    nAgents = size(obj.agents, 1);
-
-    % Network endpoints: agent 1 and the largest ODD index (nAgents when
-    % odd, nAgents-1 when even); with one or two agents only agent 1.
-    endpoints = unique([1, 2 * ceil(nAgents / 2) - 1]);
-    isEndpoint = false(1, nAgents);
-    isEndpoint(endpoints) = true;
-
-    % Trivial network: every node is an endpoint, no bulk routing needed
-    % (the lesser-neighbor constraint topology is left untouched)
-    if all(isEndpoint)
-        obj.routingAdjacencyMatrix = false(nAgents);
-        obj.routingFlows = zeros(nAgents);
-        return;
-    end
-
-    positions = zeros(nAgents, 3);
-    for kk = 1:nAgents
-        positions(kk, :) = obj.agents{kk}.pos;
-    end
-
-    % Candidate directed links: threshold-feasible pairs, excluding
-    % endpoint-sourced links (endpoints only sink; allowing them to
-    % transmit only enables useless cycles). Capacity is the Shannon
-    % spectral efficiency of the transmission direction.
-    src = zeros(nAgents^2, 1);
-    dst = zeros(nAgents^2, 1);
-    cap = zeros(nAgents^2, 1);
-    nE = 0;
-    for ii = 1:nAgents
-        if isEndpoint(ii)
-            continue;
-        end
-        for jj = 1:nAgents
-            if jj == ii || ~obj.adjacency(ii, jj)
-                continue;
-            end
-            nE = nE + 1;
-            src(nE) = ii;
-            dst(nE) = jj;
-            cap(nE) = log2(1 + obj.sinrLink(positions, jj, ii)); % tx ii -> rx jj
-        end
-    end
-    src = src(1:nE);
-    dst = dst(1:nE);
-    cap = cap(1:nE);
-
-    % --- LP over x = [f_1 .. f_nE, W] ------------------------------------
-    % min  W + eps*sum(f)   (eps breaks ties toward sparse, cycle-free flow)
-    % s.t. flow conservation (per non-endpoint i): out(i) - in(i) = 1
-    %      workload (per node i): sum(f_out/c) + sum(f_in/c) - W <= 0
-    %      f >= 0
-    relays = find(~isEndpoint);
-    nRelays = numel(relays);
-
-    fObj = [1e-6 * ones(nE, 1); 1];
-
-    Aeq = zeros(nRelays, nE + 1);
-    for rr = 1:nRelays
-        Aeq(rr, 1:nE) = double(src == relays(rr))' - double(dst == relays(rr))';
-    end
-    beq = ones(nRelays, 1);
-
-    Aineq = zeros(nAgents, nE + 1);
-    for ii = 1:nAgents
-        Aineq(ii, 1:nE) = (double(src == ii)' + double(dst == ii)') ./ cap';
-        Aineq(ii, nE + 1) = -1;
-    end
-    bineq = zeros(nAgents, 1);
-
-    lb = zeros(nE + 1, 1);
-    options = optimoptions("linprog", "Display", "off");
-
-    [x, ~, exitflag] = linprog(fObj, Aineq, bineq, Aeq, beq, lb, [], options);
-    if exitflag <= 0
-        error("Routing LP infeasible: feasible links cannot route all data to an endpoint");
-    end
-
-    % Prune trivial links and re-solve with them disallowed, so negligible
-    % flows cannot reappear in the final solution
-    ub = inf(nE + 1, 1);
-    ub(x(1:nE) < obj.routingFlowThreshold) = 0;
-    [x, ~, exitflag] = linprog(fObj, Aineq, bineq, Aeq, beq, lb, ub, options);
-    if exitflag <= 0
-        error("Routing LP infeasible after pruning: routingFlowThreshold eliminated links needed to route all data to an endpoint");
-    end
-
-    % Extract the directed topology from the pruned solution
-    flows = x(1:nE);
-    obj.routingFlows = zeros(nAgents);
-    obj.routingAdjacencyMatrix = false(nAgents);
-    for ee = 1:nE
-        if flows(ee) > 1e-9
-            obj.routingFlows(src(ee), dst(ee)) = flows(ee);
-            obj.routingAdjacencyMatrix(src(ee), dst(ee)) = true;
-        end
-    end
-
-    % OR the bulk links into the low-level (lesser-neighbor) constraint
-    % topology already on obj: the CBF maintains both tiers, keeping the
-    % network connected for basic traffic while holding the routed links.
-    % Both directions of every link are maintained (feasibility requires
-    % both to clear the threshold anyway).
-    obj.constraintAdjacencyMatrix = obj.constraintAdjacencyMatrix | obj.routingAdjacencyMatrix | obj.routingAdjacencyMatrix';
-end

@miSim/run.m

@@ -42,12 +42,13 @@ function [obj] = run(obj)
         % plotting purposes, the real partitioning is done by the agents)
         [obj.partitioning, obj.agents] = obj.agents{1}.partition(obj.agents, obj.domain.objective);
 
-        % Determine desired communications links: lesser-neighbor gives the
-        % low-level connectivity topology; routing mode overlays bulk links
+        % Determine desired communications links: the Capacity-Aware Lesser
+        % Sink Neighbor algorithm replaces lesser-neighbor when selected
         if ~obj.useFixedTopology
-            obj = obj.lesserNeighbor();
             if coder.target('MATLAB') && obj.useRoutingTopology
-                obj = obj.routeTopology();
+                obj = obj.lesserSinkNeighbor();
+            else
+                obj = obj.lesserNeighbor();
             end
         end
 

@miSim/validate.m

@@ -19,6 +19,12 @@ function validate(obj)
         error("Network is not connected");
     end
 
+    % In routing (CALSN) mode the MAINTAINED topology itself must form a
+    % single connected component (basin extension + stitch guard invariant)
+    if obj.useRoutingTopology && max(conncomp(graph(obj.constraintAdjacencyMatrix))) ~= 1
+        error("Maintained routing topology (A_control) is not connected");
+    end
+
     if any(obj.adjacency - obj.constraintAdjacencyMatrix < 0, "all")
         error("Eliminated network connections that were necessary");
     end

@miSim/writeInits.m

@@ -62,7 +62,7 @@ function writeInits(obj)
                     "useDoubleIntegrator", obj.useDoubleIntegrator, "dampingCoeff", obj.dampingCoeff, "useFixedTopology", obj.useFixedTopology, ...
                     "useSinrComms", obj.useSinrComms, "sinrThreshold", obj.sinrThreshold, "pathLossExponent", obj.pathLossExponent, "ambientTemp", obj.ambientTemp, ...
                     "centerFreq", obj.centerFreq, "bandwidth", obj.bandwidth, "txPower", txPower, ... % SINR comms model parameters
-                    "useRoutingTopology", obj.useRoutingTopology, "routingFlowThreshold", obj.routingFlowThreshold, ... % routing topology parameters
+                    "useRoutingTopology", obj.useRoutingTopology, ... % topology-selection algorithm (CALSN when true)
                     "tilt", tilt, "azimuth", azimuth, ... % joint sensor parameters
                     "alphaDist", alphaDist, "betaDist", betaDist, "alphaTilt", alphaTilt, "betaTilt", betaTilt, ... % sigmoid sensor parameters
                     "lossExponent", lossExponent, "P_TX", P_TX, "BW", BW, "f_c", f_c, "G_RX_dBi", G_RX_dBi, "beamwidthExponent", beamwidthExponent, ... % RF sensor parameters