From acfb2719e240c81dd347518c6f1d0f6d144ad45a Mon Sep 17 00:00:00 2001 From: Kevin D Date: Wed, 22 Jul 2026 17:37:04 -0700 Subject: [PATCH] new routing algorithm --- @miSim/initialize.m | 49 +++-- @miSim/initializeFromCsv.m | 9 +- @miSim/initializeFromInits.m | 7 +- @miSim/lesserSinkNeighbor.m | 372 ++++++++++++++++++++++++++++++++++ @miSim/miSim.m | 20 +- @miSim/plotGraph.m | 26 ++- @miSim/routeTopology.m | 146 -------------- @miSim/run.m | 9 +- @miSim/validate.m | 6 + @miSim/writeInits.m | 2 +- test/test_miSim.m | 376 +++++++++++++++++++++-------------- 11 files changed, 671 insertions(+), 351 deletions(-) create mode 100644 @miSim/lesserSinkNeighbor.m delete mode 100644 @miSim/routeTopology.m diff --git a/@miSim/initialize.m b/@miSim/initialize.m index 058a161..86f143b 100644 --- a/@miSim/initialize.m +++ b/@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 diff --git a/@miSim/initializeFromCsv.m b/@miSim/initializeFromCsv.m index 7929da4..6ae8151 100644 --- a/@miSim/initializeFromCsv.m +++ b/@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 diff --git a/@miSim/initializeFromInits.m b/@miSim/initializeFromInits.m index 41a7930..e0315e0 100644 --- a/@miSim/initializeFromInits.m +++ b/@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 diff --git a/@miSim/lesserSinkNeighbor.m b/@miSim/lesserSinkNeighbor.m new file mode 100644 index 0000000..ff201fd --- /dev/null +++ b/@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 diff --git a/@miSim/miSim.m b/@miSim/miSim.m index e9e27cd..e15e52f 100644 --- a/@miSim/miSim.m +++ b/@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); diff --git a/@miSim/plotGraph.m b/@miSim/plotGraph.m index 56c58ff..30f43ea 100644 --- a/@miSim/plotGraph.m +++ b/@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 diff --git a/@miSim/routeTopology.m b/@miSim/routeTopology.m deleted file mode 100644 index e642d3d..0000000 --- a/@miSim/routeTopology.m +++ /dev/null @@ -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 diff --git a/@miSim/run.m b/@miSim/run.m index f681930..c67b556 100644 --- a/@miSim/run.m +++ b/@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 diff --git a/@miSim/validate.m b/@miSim/validate.m index 370a838..9891446 100644 --- a/@miSim/validate.m +++ b/@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 diff --git a/@miSim/writeInits.m b/@miSim/writeInits.m index 33b6f39..8dae8b5 100644 --- a/@miSim/writeInits.m +++ b/@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 diff --git a/test/test_miSim.m b/test/test_miSim.m index 0c783ec..faecc1e 100644 --- a/test/test_miSim.m +++ b/test/test_miSim.m @@ -32,8 +32,8 @@ classdef test_miSim < matlab.unittest.TestCase % Agents initialStepSize = 0.2; % gradient ascent step size at the first iteration. Decreases linearly to 0 based on maxIter. initialMaxAngleStepSize = 0.1; % angular step size (degrees) for tilt/azimuth gradient ascent per timestep. - minAgents = 3; % Minimum number of agents to be randomly generated - maxAgents = 4; % Maximum number of agents to be randomly generated + minAgents = 7; % Minimum number of agents to be randomly generated + maxAgents = 7; % Maximum number of agents to be randomly generated useDoubleIntegrator = false; dampingCoeff = 2; agents = cell(0, 1); @@ -105,7 +105,7 @@ classdef test_miSim < matlab.unittest.TestCase methods (Test) % Test methods function miSim_initialization(tc) - % Test flag: set true to exercise the min-max-workload LP routing + % Test flag: set true to exercise the SINR-aware routing % overlay on this random scenario. Routing requires SINR comms, so % the flag switches the whole SINR stack on (agents keep their % default 0.1 W txPower); the permissive -30 dB threshold keeps @@ -248,7 +248,7 @@ classdef test_miSim < matlab.unittest.TestCase % Initialize the simulation if useRoutingTopology - tc.testClass = tc.testClass.initialize(tc.domain, tc.agents, tc.barrierGain, tc.barrierExponent, tc.minAlt, tc.timestep, tc.maxIter, tc.obstacles, tc.makePlots, tc.makeVideo, tc.useDoubleIntegrator, tc.dampingCoeff, tc.useFixedTopology, tc.optimizeSensorPointing, true, -30.0, 2.0, 290, 2.4e9, 20e6, true, 0.02); + tc.testClass = tc.testClass.initialize(tc.domain, tc.agents, tc.barrierGain, tc.barrierExponent, tc.minAlt, tc.timestep, tc.maxIter, tc.obstacles, tc.makePlots, tc.makeVideo, tc.useDoubleIntegrator, tc.dampingCoeff, tc.useFixedTopology, tc.optimizeSensorPointing, true, -30.0, 2.0, 290, 2.4e9, 20e6, true); else tc.testClass = tc.testClass.initialize(tc.domain, tc.agents, tc.barrierGain, tc.barrierExponent, tc.minAlt, tc.timestep, tc.maxIter, tc.obstacles, tc.makePlots, tc.makeVideo, tc.useDoubleIntegrator, tc.dampingCoeff, tc.useFixedTopology, tc.optimizeSensorPointing); end @@ -402,7 +402,7 @@ classdef test_miSim < matlab.unittest.TestCase tc.testClass = tc.testClass.run(); end function miSim_run(tc) - % Test flag: set true to exercise the min-max-workload LP routing + % Test flag: set true to exercise the min-max-workload routing % overlay on this random scenario. Routing requires SINR comms, so % the flag switches the whole SINR stack on (agents keep their % default 0.1 W txPower); the permissive -30 dB threshold keeps @@ -546,7 +546,7 @@ classdef test_miSim < matlab.unittest.TestCase % Initialize the simulation if useRoutingTopology - tc.testClass = tc.testClass.initialize(tc.domain, tc.agents, tc.barrierGain, tc.barrierExponent, tc.minAlt, tc.timestep, tc.maxIter, tc.obstacles, tc.makePlots, tc.makeVideo, tc.useDoubleIntegrator, tc.dampingCoeff, tc.useFixedTopology, tc.optimizeSensorPointing, true, -30.0, 2.0, 290, 2.4e9, 20e6, true, 0.02); + tc.testClass = tc.testClass.initialize(tc.domain, tc.agents, tc.barrierGain, tc.barrierExponent, tc.minAlt, tc.timestep, tc.maxIter, tc.obstacles, tc.makePlots, tc.makeVideo, tc.useDoubleIntegrator, tc.dampingCoeff, tc.useFixedTopology, tc.optimizeSensorPointing, true, -30.0, 2.0, 290, 2.4e9, 20e6, true); else tc.testClass = tc.testClass.initialize(tc.domain, tc.agents, tc.barrierGain, tc.barrierExponent, tc.minAlt, tc.timestep, tc.maxIter, tc.obstacles, tc.makePlots, tc.makeVideo, tc.useDoubleIntegrator, tc.dampingCoeff, tc.useFixedTopology, tc.optimizeSensorPointing); end @@ -981,129 +981,52 @@ classdef test_miSim < matlab.unittest.TestCase % Run the simulation tc.testClass = tc.testClass.run(); end - function test_routing_topology_3_agents(tc) - % Three near-collinear agents with the routing-LP topology - % selector. Endpoints are agents 1 and 3 (largest odd index); - % agent 2 is the only relay, so only its unit of data transits - % the network. Agent 2 sits slightly closer to agent 1, giving - % the 2->1 link strictly higher capacity than 2->3. The min-max - % optimum is then unique: node 2's transmit workload - % f21/c21 + f23/c23 is always the system max (endpoint receive - % time is a subset of it), so ALL flow goes over the better - % link. No split occurs because the receive side never binds - % with a single relay. - tc.minDimension = 10; % domain size + function test_routing_topology_convergence_demo(tc) + % Demonstration of the Capacity-Aware Lesser Sink Neighbor + % (CALSN) topology + routing algorithm: four agents with + % sigmoid sensors in a cube domain converge on a single + % sensing objective while the communication links evolve. + % + % Agents start scattered in a quadrilateral near one corner; + % the objective sits in the opposite corner. Base stations are + % agents 1 and 4. As the fleet migrates and contracts around + % the objective, pairwise SINRs — and with them the Shannon + % capacities, base-station potentials, parent choices, and + % capacity-proportional flow splits — change every step, so + % the network-graph tile shows the maintained (dotted) and + % flow-carrying (solid, arrowed, flow-labelled) links + % re-routing dynamically. validate() enforces the connected- + % topology invariant on every step of the run. + tc.minDimension = 10; % cube domain edge length tc.domain = tc.domain.initialize([zeros(1, 3);tc.minDimension* ones(1, 3)], REGION_TYPE.DOMAIN, "Domain"); - % make basic sensing objective - tc.domain.objective = tc.domain.objective.initialize(objectiveFunctionWrapper([2, 8; 8, 8]), tc.domain, tc.discretizationStep, tc.protectedRange); + % Single sensing objective in the far corner + tc.domain.objective = tc.domain.objective.initialize(objectiveFunctionWrapper([7, 7]), tc.domain, tc.discretizationStep, tc.protectedRange); - % Initialize agent collision geometry. Asymmetric spacing: agent 2 - % is 1.8 m from agent 1 and 2.2 m from agent 3 (both above the 1 m - % path-loss floor). The 1-3 link (4 m, with agent 2 interfering - % close-in) stays below the threshold in both directions. - tc.agents = {agent; agent; agent;}; - tc.collisionRanges = .25 * ones(size(tc.agents)); - dLeft = [1.8, 0, 0]; - dRight = [2.2, 0, 0]; - geometry1 = spherical; - geometry2 = geometry1; - geometry3 = geometry1; - geometry1 = geometry1.initialize(tc.domain.center - dLeft, tc.collisionRanges(1), REGION_TYPE.COLLISION); - geometry2 = geometry2.initialize(tc.domain.center, tc.collisionRanges(2), REGION_TYPE.COLLISION); - geometry3 = geometry3.initialize(tc.domain.center + dRight, tc.collisionRanges(3), REGION_TYPE.COLLISION); - - % Initialize agent sensor model - tc.sensor = sigmoidSensor; - tc.sensor = tc.sensor.initialize(tc.minDimension / 2, 3, 15, 3); - - % Initialize obstacles - tc.obstacles = {}; - - % Initialize agents - tc.maxIter = 50; - tc.commsRanges = 4 * ones(size(tc.agents)); % unused in SINR mode - tc.agents{1} = tc.agents{1}.initialize(tc.domain.center - dLeft, geometry1, tc.sensor, tc.commsRanges(1), tc.maxIter, tc.initialStepSize, tc.initialMaxAngleStepSize); - tc.agents{2} = tc.agents{2}.initialize(tc.domain.center, geometry2, tc.sensor, tc.commsRanges(2), tc.maxIter, tc.initialStepSize, tc.initialMaxAngleStepSize); - tc.agents{3} = tc.agents{3}.initialize(tc.domain.center + dRight, geometry3, tc.sensor, tc.commsRanges(3), tc.maxIter, tc.initialStepSize, tc.initialMaxAngleStepSize); - - % SINR communications model parameters (same threshold logic as - % the lesser-neighbor 3-agent test: a couple dB above the - % rejected end-to-end link, below the chain-link binding) - useSinrComms = true; - txPower = [0.1, 0.1, 0.1]; % homogeneous transmit power (W) - pathLossExponent = 2.0; - ambientTemp = 290; % K - centerFreq = 2.4e9; % Hz - bandwidth = 20e6; % Hz - endBindingDb = 10 * log10(0.5^pathLossExponent); % ~ -6.02 dB - sinrThreshold = endBindingDb + 2.0; % ~ -4.02 dB - useRoutingTopology = true; - routingFlowThreshold = 0.02; - - tc.agents{1}.txPower = txPower(1); - tc.agents{2}.txPower = txPower(2); - tc.agents{3}.txPower = txPower(3); - - % Initialize the simulation with the routing-LP topology selector - tc.testClass = tc.testClass.initialize(tc.domain, tc.agents, tc.barrierGain, tc.barrierExponent, tc.minAlt, tc.timestep, tc.maxIter, tc.obstacles, tc.makePlots, tc.makeVideo, tc.useDoubleIntegrator, tc.dampingCoeff, tc.useFixedTopology, tc.optimizeSensorPointing, useSinrComms, sinrThreshold, pathLossExponent, ambientTemp, centerFreq, bandwidth, useRoutingTopology, routingFlowThreshold); - - % Feasible (threshold) links form the 1-2-3 line; the routing LP - % then puts all of relay 2's data on the higher-capacity 2->1 - % link and drops 2->3 entirely - tc.assertEqual(tc.testClass.routingAdjacencyMatrix, logical( ... - [ 0, 0, 0; ... - 1, 0, 0; ... - 0, 0, 0;])); - tc.assertEqual(tc.testClass.routingFlows(2, 1), 1.0, "AbsTol", 1e-6); - - % Endpoints never transmit; total flow equals the relay's unit - tc.assertEqual(sum(tc.testClass.routingFlows, "all"), 1.0, "AbsTol", 1e-6); - - % The CBF maintains the low-level (lesser-neighbor) connectivity - % line PLUS the routed bulk link (2->1, already part of the line): - % the network must stay connected for basic traffic even though - % bulk data only flows over 2->1 - tc.assertEqual(tc.testClass.constraintAdjacencyMatrix, logical( ... - [ 1, 1, 0; ... - 1, 1, 1; ... - 0, 1, 1;])); - - % Run the simulation - tc.testClass = tc.testClass.run(); - end - function test_routing_topology_4_agents(tc) - % Four collinear agents, 2 m spacing. Endpoints are agents 1 and - % 3; relays are agents 2 and 4. Threshold-feasible links form the - % chain 1-2-3-4 (longer links are interference-limited below the - % threshold). Relay 4's only outlet is endpoint 3; relay 2 can - % reach both endpoints but 2->1 has strictly higher capacity than - % 2->3 (endpoint 3 is also loaded by relay 4's traffic), so the - % unique optimum is 2->1 and 4->3. - tc.minDimension = 10; % domain size - tc.domain = tc.domain.initialize([zeros(1, 3);tc.minDimension* ones(1, 3)], REGION_TYPE.DOMAIN, "Domain"); - - % make basic sensing objective - tc.domain.objective = tc.domain.objective.initialize(objectiveFunctionWrapper([2, 8; 8, 8]), tc.domain, tc.discretizationStep, tc.protectedRange); - - % Initialize agent collision geometry: collinear at x = 2, 4, 6, 8 + % Four agents scattered (non-collinear) near the (low-x, low-y) + % corner, co-altitude at z = 2. Base 1 trails the fleet and + % base 4 leads it, with the relays staggered between: relay 3 + % starts too far from base 1 for a direct link to be worth + % using, so its initial basin-1 route is MULTI-HOP via relay 2 + % and re-parents to the base directly as the fleet contracts. tc.agents = {agent; agent; agent; agent;}; tc.collisionRanges = .25 * ones(size(tc.agents)); - d = [2.0, 0, 0]; - positions = [tc.domain.center - 1.5 * d; ... - tc.domain.center - 0.5 * d; ... - tc.domain.center + 0.5 * d; ... - tc.domain.center + 1.5 * d;]; + positions = [1.0, 2.0, 2; ... + 3.0, 3.5, 2; ... + 5.5, 3.0, 2; ... + 5.0, 5.5, 2;]; - % Initialize agent sensor model + % Initialize agent sensor model (sigmoid) tc.sensor = sigmoidSensor; tc.sensor = tc.sensor.initialize(tc.minDimension / 2, 3, 15, 3); % Initialize obstacles tc.obstacles = {}; - % Initialize agents - tc.maxIter = 50; + % Initialize agents: 50 steps is long enough to watch the + % fleet converge on the objective and settle (the gradient + % step decays linearly to zero over maxIter) + tc.maxIter = 80; tc.commsRanges = 4 * ones(size(tc.agents)); % unused in SINR mode for aa = 1:4 geometry = spherical; @@ -1112,48 +1035,195 @@ classdef test_miSim < matlab.unittest.TestCase tc.agents{aa}.txPower = 0.1; end - % SINR communications model parameters - useSinrComms = true; - pathLossExponent = 2.0; - ambientTemp = 290; % K - centerFreq = 2.4e9; % Hz - bandwidth = 20e6; % Hz - endBindingDb = 10 * log10(0.5^pathLossExponent); % ~ -6.02 dB - sinrThreshold = endBindingDb + 2.0; % ~ -4.02 dB - useRoutingTopology = true; - routingFlowThreshold = 0.02; + % SINR communications model. The threshold must sit >= 3 dB + % below the worst link CALSN may select, because the CBF holds + % every maintained link at threshold + 3 dB in BOTH directions: + % a link selected inside that 3 dB band is born with a violated + % barrier and can deadlock the QP (agents hold position, the + % geometry never changes, and the same link is re-selected + % forever). The staggered start's weakest candidate link (the + % long 1-4 base-to-base diagonal) sits at ~ -9 dB and the + % weakest link CALSN actually selects at ~ -4 dB, so -14 dB + % leaves the -11 dB CBF floor a comfortable margin. + useSinrComms = true; + pathLossExponent = 2.0; + ambientTemp = 290; % K + centerFreq = 2.4e9; % Hz + bandwidth = 20e6; % Hz + sinrThreshold = -14.0; % dB + useRoutingTopology = true; - % Initialize the simulation with the routing-LP topology selector - tc.testClass = tc.testClass.initialize(tc.domain, tc.agents, tc.barrierGain, tc.barrierExponent, tc.minAlt, tc.timestep, tc.maxIter, tc.obstacles, tc.makePlots, tc.makeVideo, tc.useDoubleIntegrator, tc.dampingCoeff, tc.useFixedTopology, tc.optimizeSensorPointing, useSinrComms, sinrThreshold, pathLossExponent, ambientTemp, centerFreq, bandwidth, useRoutingTopology, routingFlowThreshold); + % Initialize the simulation with the CALSN topology selector + tc.testClass = tc.testClass.initialize(tc.domain, tc.agents, tc.barrierGain, tc.barrierExponent, tc.minAlt, tc.timestep, tc.maxIter, tc.obstacles, tc.makePlots, tc.makeVideo, tc.useDoubleIntegrator, tc.dampingCoeff, tc.useFixedTopology, tc.optimizeSensorPointing, useSinrComms, sinrThreshold, pathLossExponent, ambientTemp, centerFreq, bandwidth, useRoutingTopology); - % Unique optimum: 2->1 and 4->3 - expectedDirected = false(4); - expectedDirected(2, 1) = true; - expectedDirected(4, 3) = true; - tc.assertEqual(tc.testClass.routingAdjacencyMatrix, expectedDirected); - - % The CBF maintains the full lesser-neighbor chain (low-level - % connectivity) with the bulk links overlaid on it - tc.assertEqual(tc.testClass.constraintAdjacencyMatrix, logical( ... - [ 1, 1, 0, 0; ... - 1, 1, 1, 0; ... - 0, 1, 1, 1; ... - 0, 0, 1, 1;])); - - % Flow conservation at every relay: out - in = 1 unit produced + % Initial routing invariants: bases at zero potential, relays + % finite; maintained topology one connected component; both + % generated units delivered to the bases with nothing unrouted + phi = tc.testClass.routingPotentials; + tc.assertEqual(phi([1, 4]), [0; 0]); + tc.assertTrue(all(isfinite(phi([2, 3]))) && all(phi([2, 3]) > 0)); + tc.assertEqual(max(conncomp(graph(tc.testClass.constraintAdjacencyMatrix))), 1); F = tc.testClass.routingFlows; - endpoints = [1, 3]; - for aa = 1:4 - if any(aa == endpoints) - continue; + for aa = [2, 3] + tc.assertEqual(sum(F(aa, :)) - sum(F(:, aa)), 1.0, "AbsTol", 1e-9); + end + bases = [1, 4]; + tc.assertEqual(sum(F(:, bases), "all") - sum(F(bases, :), "all"), 2.0, "AbsTol", 1e-9); + tc.assertTrue(all(tc.testClass.routingUnrouted == 0)); + + % Run the simulation (per-step CALSN re-solve; validate() + % checks feasibility- and maintained-graph connectivity at + % every iteration) + tc.testClass = tc.testClass.run(); + + % Convergence: total sensing performance improved over the run + tc.assertGreaterThan(tc.testClass.performance(end), tc.testClass.performance(2)); + + % Dynamic evolution: the maintained topology changed at least + % once relative to the initial configuration + hist = tc.testClass.constraintAdjacencyHist; + changed = false; + for kk = 2:size(hist, 3) + if any(hist(:, :, kk) ~= hist(:, :, 1), "all") + changed = true; + break; end - tc.assertEqual(sum(F(aa, :)) - sum(F(:, aa)), 1.0, "AbsTol", 1e-6); + end + tc.assertTrue(changed, "Maintained topology never changed during the run"); + + % Final routing state is still a valid delivery tree: one + % connected component, conservation at the generating UAVs, + % both units delivered, nothing unrouted + tc.assertEqual(max(conncomp(graph(tc.testClass.constraintAdjacencyMatrix))), 1); + F = tc.testClass.routingFlows; + for aa = [2, 3] + tc.assertEqual(sum(F(aa, :)) - sum(F(:, aa)), 1.0, "AbsTol", 1e-9); + end + tc.assertEqual(sum(F(:, bases), "all") - sum(F(bases, :), "all"), 2.0, "AbsTol", 1e-9); + tc.assertTrue(all(tc.testClass.routingUnrouted == 0)); + end + function test_routing_topology_convergence_demo_rfSensor(tc) + % Demonstration of the Capacity-Aware Lesser Sink Neighbor + % (CALSN) topology + routing algorithm: four agents with + % sigmoid sensors in a cube domain converge on a single + % sensing objective while the communication links evolve. + % + % Agents start scattered in a quadrilateral near one corner; + % the objective sits in the opposite corner. Base stations are + % agents 1 and 4. As the fleet migrates and contracts around + % the objective, pairwise SINRs — and with them the Shannon + % capacities, base-station potentials, parent choices, and + % capacity-proportional flow splits — change every step, so + % the network-graph tile shows the maintained (dotted) and + % flow-carrying (solid, arrowed, flow-labelled) links + % re-routing dynamically. validate() enforces the connected- + % topology invariant on every step of the run. + tc.minDimension = 10; % cube domain edge length + tc.domain = tc.domain.initialize([zeros(1, 3);tc.minDimension* ones(1, 3)], REGION_TYPE.DOMAIN, "Domain"); + + % Single sensing objective in the far corner + tc.domain.objective = tc.domain.objective.initialize(objectiveFunctionWrapper([7, 7]), tc.domain, tc.discretizationStep, tc.protectedRange); + + % Four agents scattered (non-collinear) near the (low-x, low-y) + % corner, co-altitude at z = 2. Base 1 trails the fleet and + % base 4 leads it, with the relays staggered between: relay 3 + % starts too far from base 1 for a direct link to be worth + % using, so its initial basin-1 route is MULTI-HOP via relay 2 + % and re-parents to the base directly as the fleet contracts. + tc.agents = {agent; agent; agent; agent;}; + tc.collisionRanges = .25 * ones(size(tc.agents)); + positions = [1.0, 2.0, 2; ... + 3.0, 3.5, 2; ... + 5.5, 3.0, 2; ... + 5.0, 5.5, 2;]; + + % Initialize agent sensor model (sigmoid) + tc.sensor = rfSensor; + tilt = 0; azimuth = 0; + tc.sensor = tc.sensor.initialize(tc.P_TX * 1 + rand * 4, tc.T_0, tc.BW, tc.f_c, tc.G_RX_dBi, tc.beamwidthExponent + randi(100), tilt, azimuth, tc.lossExponent); + + % Initialize obstacles + tc.obstacles = {}; + + % Initialize agents: 50 steps is long enough to watch the + % fleet converge on the objective and settle (the gradient + % step decays linearly to zero over maxIter) + tc.maxIter = 80; + tc.commsRanges = 4 * ones(size(tc.agents)); % unused in SINR mode + for aa = 1:4 + geometry = spherical; + geometry = geometry.initialize(positions(aa, :), tc.collisionRanges(aa), REGION_TYPE.COLLISION); + tc.agents{aa} = tc.agents{aa}.initialize(positions(aa, :), geometry, tc.sensor, tc.commsRanges(aa), tc.maxIter, tc.initialStepSize, tc.initialMaxAngleStepSize); + tc.agents{aa}.txPower = 0.1; end - % Every surviving link carries at least the pruning threshold, - % and total delivered data equals the number of relays - tc.assertTrue(all(F(F > 0) >= tc.testClass.routingFlowThreshold)); - tc.assertEqual(sum(F(:, endpoints), "all") - sum(F(endpoints, :), "all"), 2.0, "AbsTol", 1e-6); + % SINR communications model. The threshold must sit >= 3 dB + % below the worst link CALSN may select, because the CBF holds + % every maintained link at threshold + 3 dB in BOTH directions: + % a link selected inside that 3 dB band is born with a violated + % barrier and can deadlock the QP (agents hold position, the + % geometry never changes, and the same link is re-selected + % forever). The staggered start's weakest candidate link (the + % long 1-4 base-to-base diagonal) sits at ~ -9 dB and the + % weakest link CALSN actually selects at ~ -4 dB, so -14 dB + % leaves the -11 dB CBF floor a comfortable margin. + useSinrComms = true; + pathLossExponent = 2.0; + ambientTemp = 290; % K + centerFreq = 2.4e9; % Hz + bandwidth = 20e6; % Hz + sinrThreshold = -14.0; % dB + useRoutingTopology = true; + + % Initialize the simulation with the CALSN topology selector + tc.optimizeSensorPointing = true; + tc.testClass = tc.testClass.initialize(tc.domain, tc.agents, tc.barrierGain, tc.barrierExponent, tc.minAlt, tc.timestep, tc.maxIter, tc.obstacles, tc.makePlots, tc.makeVideo, tc.useDoubleIntegrator, tc.dampingCoeff, tc.useFixedTopology, tc.optimizeSensorPointing, useSinrComms, sinrThreshold, pathLossExponent, ambientTemp, centerFreq, bandwidth, useRoutingTopology); + + % Initial routing invariants: bases at zero potential, relays + % finite; maintained topology one connected component; both + % generated units delivered to the bases with nothing unrouted + phi = tc.testClass.routingPotentials; + tc.assertEqual(phi([1, 4]), [0; 0]); + tc.assertTrue(all(isfinite(phi([2, 3]))) && all(phi([2, 3]) > 0)); + tc.assertEqual(max(conncomp(graph(tc.testClass.constraintAdjacencyMatrix))), 1); + F = tc.testClass.routingFlows; + for aa = [2, 3] + tc.assertEqual(sum(F(aa, :)) - sum(F(:, aa)), 1.0, "AbsTol", 1e-9); + end + bases = [1, 4]; + tc.assertEqual(sum(F(:, bases), "all") - sum(F(bases, :), "all"), 2.0, "AbsTol", 1e-9); + tc.assertTrue(all(tc.testClass.routingUnrouted == 0)); + + % Run the simulation (per-step CALSN re-solve; validate() + % checks feasibility- and maintained-graph connectivity at + % every iteration) + tc.testClass = tc.testClass.run(); + + % Convergence: total sensing performance improved over the run + tc.assertGreaterThan(tc.testClass.performance(end), tc.testClass.performance(2)); + + % Dynamic evolution: the maintained topology changed at least + % once relative to the initial configuration + hist = tc.testClass.constraintAdjacencyHist; + changed = false; + for kk = 2:size(hist, 3) + if any(hist(:, :, kk) ~= hist(:, :, 1), "all") + changed = true; + break; + end + end + tc.assertTrue(changed, "Maintained topology never changed during the run"); + + % Final routing state is still a valid delivery tree: one + % connected component, conservation at the generating UAVs, + % both units delivered, nothing unrouted + tc.assertEqual(max(conncomp(graph(tc.testClass.constraintAdjacencyMatrix))), 1); + F = tc.testClass.routingFlows; + for aa = [2, 3] + tc.assertEqual(sum(F(aa, :)) - sum(F(:, aa)), 1.0, "AbsTol", 1e-9); + end + tc.assertEqual(sum(F(:, bases), "all") - sum(F(bases, :), "all"), 2.0, "AbsTol", 1e-9); + tc.assertTrue(all(tc.testClass.routingUnrouted == 0)); end function test_communications_constraint_fixed_radius(tc) % No obstacles