From 37598479cc7b0f35952baf4a914cef5f5828353a Mon Sep 17 00:00:00 2001 From: Kevin D Date: Sun, 12 Jul 2026 22:50:14 -0700 Subject: [PATCH] better routing solution --- @miSim/initialize.m | 18 +++- @miSim/initializeFromCsv.m | 15 ++- @miSim/initializeFromInits.m | 13 ++- @miSim/miSim.m | 11 +- @miSim/plotGraph.m | 116 ++++++++++---------- @miSim/routeTopology.m | 146 +++++++++++++++++++++++++ @miSim/run.m | 6 +- @miSim/validate.m | 9 ++ @miSim/writeInits.m | 1 + test/test_miSim.m | 200 ++++++++++++++++++++++++++++++++++- 10 files changed, 469 insertions(+), 66 deletions(-) create mode 100644 @miSim/routeTopology.m diff --git a/@miSim/initialize.m b/@miSim/initialize.m index e51cb75..058a161 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) +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) arguments (Input) obj (1, 1) {mustBeA(obj, "miSim")}; domain (1, 1) {mustBeGeometry}; @@ -21,6 +21,8 @@ function [obj] = initialize(obj, domain, agents, barrierGain, barrierExponent, m ambientTemp (1, 1) double = 290.0; 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")}; @@ -112,12 +114,24 @@ function [obj] = initialize(obj, domain, agents, barrierGain, barrierExponent, m obj.centerFreq = centerFreq; obj.bandwidth = bandwidth; - % Compute adjacency matrix and network topology + % Set topology-selection algorithm and routing parameters + 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. obj = obj.updateAdjacency(); if obj.useFixedTopology obj.constraintAdjacencyMatrix = obj.adjacency; else obj = obj.lesserNeighbor(); + if coder.target('MATLAB') && obj.useRoutingTopology + obj = obj.routeTopology(); + end end % Set up times to iterate over diff --git a/@miSim/initializeFromCsv.m b/@miSim/initializeFromCsv.m index 310ddf7..7929da4 100644 --- a/@miSim/initializeFromCsv.m +++ b/@miSim/initializeFromCsv.m @@ -147,6 +147,18 @@ else BANDWIDTH = 20e6; end +% Routing topology selection (optional columns; 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; dom = dom.initialize([DOMAIN_MIN; DOMAIN_MAX], REGION_TYPE.DOMAIN, "Guidance Domain"); @@ -195,6 +207,7 @@ end 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_SINR_COMMS, SINR_THRESHOLD, PATH_LOSS_EXPONENT, AMBIENT_TEMP, CENTER_FREQ, BANDWIDTH, ... + USE_ROUTING_TOPOLOGY, ROUTING_FLOW_THRESHOLD); end diff --git a/@miSim/initializeFromInits.m b/@miSim/initializeFromInits.m index ff3300d..41a7930 100644 --- a/@miSim/initializeFromInits.m +++ b/@miSim/initializeFromInits.m @@ -111,11 +111,22 @@ if isfield(inits, 'bandwidth') else bandwidth = 20e6; end +if isfield(inits, 'useRoutingTopology') + useRoutingTopology = logical(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); + useSinrComms, sinrThreshold, pathLossExponent, ambientTemp, centerFreq, bandwidth, ... + useRoutingTopology, routingFlowThreshold); end diff --git a/@miSim/miSim.m b/@miSim/miSim.m index 3002e9f..e9e27cd 100644 --- a/@miSim/miSim.m +++ b/@miSim/miSim.m @@ -30,6 +30,14 @@ classdef miSim ambientTemp = 290.0; % ambient temperature (K) for thermal noise N = k_B*T*B 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. + 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) artifactName = ""; f; % main plotting tiled layout figure fPerf; % performance plot figure @@ -79,12 +87,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); + [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] = initializeFromCsv(obj, csvPath); [obj] = initializeFromInits(obj, initsPath); [obj] = plotFromSimHist(obj, initsPath, histPath); [obj] = run(obj); [obj] = lesserNeighbor(obj); + [obj] = routeTopology(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 d303833..56c58ff 100644 --- a/@miSim/plotGraph.m +++ b/@miSim/plotGraph.m @@ -6,7 +6,12 @@ function obj = plotGraph(obj) obj (1, 1) {mustBeA(obj, "miSim")}; end - nAgents = size(obj.agents, 1); + % 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. + G = graph(obj.constraintAdjacencyMatrix, "omitselfloops"); % Resolve the target axes for the network graph tile if isnan(obj.networkGraphIndex) @@ -15,72 +20,67 @@ function obj = plotGraph(obj) ax = obj.f.Children(1).Children(obj.networkGraphIndex(1)); end - % Build the graph to plot. The fixed-radius model is undirected. In SINR - % mode the link set (lesser-neighbor topology) is kept exactly as-is, but - % each connection is given an orientation: of the two directional SINRs on - % a link, the larger one wins, and the edge is directed from that - % direction's transmitter to its receiver (arrow = transmitter -> receiver). + % Plot graph object(s) + 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 + % 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); + o = [o; plot(ax, D, "XData", o(1).XData, "YData", o(1).YData, ... + "LineStyle", "-", "EdgeColor", "g", "LineWidth", 2, ... + "ArrowSize", 12, "Marker", "none", "NodeLabel", {})]; + end + hold(ax, "off"); + nBase = size(o, 1); + if ~isnan(obj.networkGraphIndex) && size(obj.networkGraphIndex, 2) > 1 + for ii = 2:size(obj.networkGraphIndex, 2) + o = [o; copyobj(o(1:nBase), obj.f.Children(1).Children(obj.networkGraphIndex(ii)))]; + end + end + + % In SINR comms mode, label each edge with the connection's SINR and report + % the connectivity threshold in the tile title. The fixed-radius model is + % left unchanged (no edge labels, title stays as set in firstPlotSetup). if obj.useSinrComms + endNodes = G.Edges.EndNodes; + nAgents = size(obj.agents, 1); positions = zeros(nAgents, 3); for kk = 1:nAgents positions(kk, :) = obj.agents{kk}.pos; end - D = false(nAgents); % directed adjacency: D(tx, rx) = true - for aa = 1:(nAgents - 1) - for bb = (aa + 1):nAgents - if obj.constraintAdjacencyMatrix(aa, bb) - sinrARx = obj.sinrLink(positions, aa, bb); % aa receives bb - sinrBRx = obj.sinrLink(positions, bb, aa); % bb receives aa - if sinrARx >= sinrBRx - D(bb, aa) = true; % better as bb -> aa (tx bb, rx aa) - else - D(aa, bb) = true; % better as aa -> bb (tx aa, rx bb) - end - end - end - end - G = digraph(D); - else - G = graph(obj.constraintAdjacencyMatrix, "omitselfloops"); - end - - % Plot graph object (digraph renders arrowheads automatically) - hold(ax, "on"); - o = plot(ax, G, "LineStyle", "--", "EdgeColor", "g", "NodeColor", "k", "LineWidth", 2); - hold(ax, "off"); - if ~isnan(obj.networkGraphIndex) && size(obj.networkGraphIndex, 2) > 1 - for ii = 2:size(obj.networkGraphIndex, 2) - o = [o; copyobj(o(1), obj.f.Children(1).Children(obj.networkGraphIndex(ii)))]; - end - end - - % Mark ground-station-connected nodes with an "x" (others keep the default - % dot). Ground stations attach to the two "router" agents: agent 1 and the - % agent with the largest ODD index (that is nAgents when odd, else - % nAgents-1). With one or two agents only agent 1 qualifies. - routers = unique([1, 2 * ceil(nAgents / 2) - 1]); - for oo = 1:numel(o) - highlight(o(oo), routers, "Marker", "x", "MarkerSize", 8); - end - - % In SINR comms mode, label each directed edge with the better SINR that set - % its orientation and report the threshold in the tile title. The - % fixed-radius model is left unchanged (no edge labels, title from setup). - if obj.useSinrComms - endNodes = G.Edges.EndNodes; % each row [tx rx] for a directed edge - % Label with the SINR of the chosen (better) direction: the receiver is - % the edge target, the transmitter the edge source. Bare numbers (no - % units) keep the small tile uncluttered; GraphPlot draws labels - % horizontally, left-to-right. + % Labels are bare numbers (no units) so they stay compact; GraphPlot + % renders them horizontally, left-to-right. Dotted (low-level) 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. edgeLabels = strings(size(endNodes, 1), 1); for ee = 1:size(endNodes, 1) - sinrLin = obj.sinrLink(positions, endNodes(ee, 2), endNodes(ee, 1)); % rx, tx - edgeLabels(ee) = sprintf("%.1f", 10 * log10(sinrLin)); + lo = endNodes(ee, 1); + hi = endNodes(ee, 2); + if obj.useRoutingTopology && (obj.routingAdjacencyMatrix(lo, hi) || obj.routingAdjacencyMatrix(hi, lo)) + edgeLabels(ee) = ""; % solid overlay carries this link's label + else + sinrLoRx = obj.sinrLink(positions, lo, hi); % receiver lo (lower index) + sinrHiRx = obj.sinrLink(positions, hi, lo); % receiver hi (higher index) + edgeLabels(ee) = sprintf("%.1f", 10 * log10(min(sinrLoRx, sinrHiRx))); + end end o(1).EdgeLabel = edgeLabels; - % Threshold in the existing single-line title (not a subtitle) so the - % small tile isn't compressed, and disable the hover toolbar so the - % R2026 "..." menu stops covering the readout. + if obj.useRoutingTopology + endNodesD = D.Edges.EndNodes; % each row [tx rx] + bulkLabels = 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)); + end + o(2).EdgeLabel = bulkLabels; + 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 + % disabled so the R2026 "..." menu stops covering the readout. title(ax, sprintf("Network Graph (\\geq %.1f dB)", obj.sinrThreshold)); ax.Toolbar.Visible = "off"; end diff --git a/@miSim/routeTopology.m b/@miSim/routeTopology.m new file mode 100644 index 0000000..e642d3d --- /dev/null +++ b/@miSim/routeTopology.m @@ -0,0 +1,146 @@ +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 aa0ed7e..f681930 100644 --- a/@miSim/run.m +++ b/@miSim/run.m @@ -42,9 +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 + % Determine desired communications links: lesser-neighbor gives the + % low-level connectivity topology; routing mode overlays bulk links if ~obj.useFixedTopology obj = obj.lesserNeighbor(); + if coder.target('MATLAB') && obj.useRoutingTopology + obj = obj.routeTopology(); + end end % Log constraint adjacency for this timestep diff --git a/@miSim/validate.m b/@miSim/validate.m index f57249e..370a838 100644 --- a/@miSim/validate.m +++ b/@miSim/validate.m @@ -5,6 +5,15 @@ function validate(obj) arguments (Output) end + %% Mode Validators + if obj.useRoutingTopology && ~obj.useSinrComms + error("Routing topology requires SINR comms mode (useSinrComms)"); + end + + if obj.useRoutingTopology && obj.useFixedTopology + error("useRoutingTopology and useFixedTopology are mutually exclusive topology selectors"); + end + %% Communications Network Validators if max(conncomp(graph(obj.adjacency))) ~= 1 error("Network is not connected"); diff --git a/@miSim/writeInits.m b/@miSim/writeInits.m index 28a9e96..33b6f39 100644 --- a/@miSim/writeInits.m +++ b/@miSim/writeInits.m @@ -62,6 +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 "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 0bde39e..0c783ec 100644 --- a/test/test_miSim.m +++ b/test/test_miSim.m @@ -105,6 +105,13 @@ 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 + % 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 + % random geometries fully connected. + useRoutingTopology = true; + % randomly create obstacles nGeom = tc.minNumObstacles + randi(tc.maxNumObstacles - tc.minNumObstacles); tc.obstacles = cell(nGeom, 1); @@ -240,7 +247,11 @@ classdef test_miSim < matlab.unittest.TestCase end % Initialize the simulation - 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); + 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); + 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 end function miSim_run_rf_sensor(tc) % randomly create obstacles @@ -391,6 +402,13 @@ 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 + % 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 + % random geometries fully connected. + useRoutingTopology = true; + % randomly create obstacles nGeom = tc.minNumObstacles + randi(tc.maxNumObstacles - tc.minNumObstacles); tc.obstacles = cell(nGeom, 1); @@ -527,7 +545,11 @@ classdef test_miSim < matlab.unittest.TestCase end % Initialize the simulation - 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); + 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); + 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 % Write out initialization state tc.testClass.writeInits(); @@ -959,6 +981,180 @@ 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 + 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. 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 + 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;]; + + % 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 + 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 + + % 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; + + % 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); + + % 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 + F = tc.testClass.routingFlows; + endpoints = [1, 3]; + for aa = 1:4 + if any(aa == endpoints) + continue; + end + tc.assertEqual(sum(F(aa, :)) - sum(F(:, aa)), 1.0, "AbsTol", 1e-6); + 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); + end function test_communications_constraint_fixed_radius(tc) % No obstacles % Fixed two agents initial conditions