7 Commits

209 changed files with 1570 additions and 1219 deletions

40
@agent/agent.m Normal file
View File

@@ -0,0 +1,40 @@
classdef agent
properties (SetAccess = private, GetAccess = public)
% Identifiers
index = NaN;
label = "";
% Sensor
sensorModel;
sensingLength = 0.05; % length parameter used by sensing function
% Guidance
guidanceModel;
% State
lastPos = NaN(1, 3); % position from previous timestep
pos = NaN(1, 3); % current position
vel = NaN(1, 3); % current velocity
pan = NaN; % pan angle
tilt = NaN; % tilt angle
% Collision
collisionGeometry;
% FOV cone
fovGeometry;
% Communication
comRange = NaN;
% Plotting
scatterPoints;
end
methods (Access = public)
[obj] = initialize(obj, pos, vel, pan, tilt, collisionGeometry, sensorModel, guidanceModel, comRange, index, label);
[obj] = run(obj, sensingObjective, domain, partitioning);
[obj, f] = plot(obj, ind, f);
updatePlots(obj);
end
end

33
@agent/initialize.m Normal file
View File

@@ -0,0 +1,33 @@
function obj = initialize(obj, pos, vel, pan, tilt, collisionGeometry, sensorModel, guidanceModel, comRange, index, label)
arguments (Input)
obj (1, 1) {mustBeA(obj, 'agent')};
pos (1, 3) double;
vel (1, 3) double;
pan (1, 1) double;
tilt (1, 1) double;
collisionGeometry (1, 1) {mustBeGeometry};
sensorModel (1, 1) {mustBeSensor}
guidanceModel (1, 1) {mustBeA(guidanceModel, 'function_handle')};
comRange (1, 1) double = NaN;
index (1, 1) double = NaN;
label (1, 1) string = "";
end
arguments (Output)
obj (1, 1) {mustBeA(obj, 'agent')};
end
obj.pos = pos;
obj.vel = vel;
obj.pan = pan;
obj.tilt = tilt;
obj.collisionGeometry = collisionGeometry;
obj.sensorModel = sensorModel;
obj.guidanceModel = guidanceModel;
obj.comRange = comRange;
obj.index = index;
obj.label = label;
% Initialize FOV cone
obj.fovGeometry = cone;
obj.fovGeometry = obj.fovGeometry.initialize([obj.pos(1:2), 0], tan(obj.sensorModel.alphaTilt) * obj.pos(3), obj.pos(3), REGION_TYPE.FOV, sprintf("%s FOV", obj.label));
end

35
@agent/plot.m Normal file
View File

@@ -0,0 +1,35 @@
function [obj, f] = plot(obj, ind, f)
arguments (Input)
obj (1, 1) {mustBeA(obj, 'agent')};
ind (1, :) double = NaN;
f (1, 1) {mustBeA(f, 'matlab.ui.Figure')} = figure;
end
arguments (Output)
obj (1, 1) {mustBeA(obj, 'agent')};
f (1, 1) {mustBeA(f, 'matlab.ui.Figure')};
end
% Create axes if they don't already exist
f = firstPlotSetup(f);
% Plot points representing the agent position
hold(f.Children(1).Children(end), "on");
o = scatter3(f.Children(1).Children(end), obj.pos(1), obj.pos(2), obj.pos(3), 'filled', 'ko', 'SizeData', 25);
hold(f.Children(1).Children(end), "off");
% Check if this is a tiled layout figure
if strcmp(f.Children(1).Type, 'tiledlayout')
% Add to other perspectives
o = [o; copyobj(o(1), f.Children(1).Children(2))];
o = [o; copyobj(o(1), f.Children(1).Children(3))];
o = [o; copyobj(o(1), f.Children(1).Children(4))];
end
obj.scatterPoints = o;
% Plot collision geometry
[obj.collisionGeometry, f] = obj.collisionGeometry.plotWireframe(ind, f);
% Plot FOV geometry
[obj.fovGeometry, f] = obj.fovGeometry.plot(ind, f);
end

28
@agent/run.m Normal file
View File

@@ -0,0 +1,28 @@
function obj = run(obj, sensingObjective, domain, partitioning)
arguments (Input)
obj (1, 1) {mustBeA(obj, 'agent')};
sensingObjective (1, 1) {mustBeA(sensingObjective, 'sensingObjective')};
domain (1, 1) {mustBeGeometry};
partitioning (:, :) double;
end
arguments (Output)
obj (1, 1) {mustBeA(obj, 'agent')};
end
% Do sensing
[sensedValues, sensedPositions] = obj.sensorModel.sense(obj, sensingObjective, domain, partitioning);
% Determine next planned position
nextPos = obj.guidanceModel(sensedValues, sensedPositions, obj.pos);
% Move to next position
% (dynamics not modeled at this time)
obj.lastPos = obj.pos;
obj.pos = nextPos;
% Calculate movement
d = obj.pos - obj.collisionGeometry.center;
% Reinitialize collision geometry in the new position
obj.collisionGeometry = obj.collisionGeometry.initialize([obj.collisionGeometry.minCorner; obj.collisionGeometry.maxCorner] + d, obj.collisionGeometry.tag, obj.collisionGeometry.label);
end

35
@agent/updatePlots.m Normal file
View File

@@ -0,0 +1,35 @@
function updatePlots(obj)
arguments (Input)
obj (1, 1) {mustBeA(obj, 'agent')};
end
arguments (Output)
end
% Scatterplot point positions
for ii = 1:size(obj.scatterPoints, 1)
obj.scatterPoints(ii).XData = obj.pos(1);
obj.scatterPoints(ii).YData = obj.pos(2);
obj.scatterPoints(ii).ZData = obj.pos(3);
end
% Find change in agent position since last timestep
deltaPos = obj.pos - obj.lastPos;
% Collision geometry edges
for jj = 1:size(obj.collisionGeometry.lines, 2)
% Update plotting
for ii = 1:size(obj.collisionGeometry.lines(:, jj), 1)
obj.collisionGeometry.lines(ii, jj).XData = obj.collisionGeometry.lines(ii, jj).XData + deltaPos(1);
obj.collisionGeometry.lines(ii, jj).YData = obj.collisionGeometry.lines(ii, jj).YData + deltaPos(2);
obj.collisionGeometry.lines(ii, jj).ZData = obj.collisionGeometry.lines(ii, jj).ZData + deltaPos(3);
end
end
% Update FOV geometry surfaces
for jj = 1:size(obj.fovGeometry.surface, 2)
% Update each plot
obj.fovGeometry.surface(jj).XData = obj.fovGeometry.surface(jj).XData + deltaPos(1);
obj.fovGeometry.surface(jj).YData = obj.fovGeometry.surface(jj).YData + deltaPos(2);
obj.fovGeometry.surface(jj).ZData = obj.fovGeometry.surface(jj).ZData + deltaPos(3);
end
end

42
@miSim/initialize.m Normal file
View File

@@ -0,0 +1,42 @@
function [obj, f] = initialize(obj, domain, objective, agents, timestep, partitoningFreq, maxIter, obstacles)
arguments (Input)
obj (1, 1) {mustBeA(obj, 'miSim')};
domain (1, 1) {mustBeGeometry};
objective (1, 1) {mustBeA(objective, 'sensingObjective')};
agents (:, 1) cell;
timestep (:, 1) double = 0.05;
partitoningFreq (:, 1) double = 0.25
maxIter (:, 1) double = 1000;
obstacles (:, 1) cell {mustBeGeometry} = cell(0, 1);
end
arguments (Output)
obj (1, 1) {mustBeA(obj, 'miSim')};
f (1, 1) {mustBeA(f, 'matlab.ui.Figure')};
end
% Define simulation time parameters
obj.timestep = timestep;
obj.maxIter = maxIter;
% Define domain
obj.domain = domain;
obj.partitioningFreq = partitoningFreq;
% Add geometries representing obstacles within the domain
obj.obstacles = obstacles;
% Define objective
obj.objective = objective;
% Define agents
obj.agents = agents;
% Compute adjacency matrix
obj = obj.updateAdjacency();
% Create initial partitioning
obj = obj.partition();
% Set up plots showing initialized state
[obj, f] = obj.plot();
end

45
@miSim/miSim.m Normal file
View File

@@ -0,0 +1,45 @@
classdef miSim
% multiagent interconnection simulation
% Simulation parameters
properties (SetAccess = private, GetAccess = public)
timestep = NaN; % delta time interval for simulation iterations
partitioningFreq = NaN; % number of simulation timesteps at which the partitioning routine is re-run
maxIter = NaN; % maximum number of simulation iterations
domain = rectangularPrism;
objective = sensingObjective;
obstacles = cell(0, 1); % geometries that define obstacles within the domain
agents = cell(0, 1); % agents that move within the domain
adjacency = NaN; % Adjacency matrix representing communications network graph
sensorPerformanceMinimum = 1e-6; % minimum sensor performance to allow assignment of a point in the domain to a partition
partitioning = NaN;
end
properties (Access = private)
% Plot objects
connectionsPlot; % objects for lines connecting agents in spatial plots
graphPlot; % objects for abstract network graph plot
partitionPlot; % objects for partition plot
% Indicies for various plot types in the main tiled layout figure
spatialPlotIndices = [6, 4, 3, 2];
objectivePlotIndices = [6, 4];
networkGraphIndex = 5;
partitionGraphIndex = 1;
end
methods (Access = public)
[obj, f] = initialize(obj, domain, objective, agents, timestep, partitoningFreq, maxIter, obstacles);
[obj, f] = run(obj, f);
[obj] = partition(obj);
[obj] = updateAdjacency(obj);
[obj, f] = plot(obj);
[obj, f] = plotConnections(obj, ind, f);
[obj, f] = plotPartitions(obj, ind, f);
[obj, f] = plotGraph(obj, ind, f);
[obj, f] = updatePlots(obj, f, updatePartitions);
end
methods (Access = private)
[v] = setupVideoWriter(obj);
end
end

27
@miSim/partition.m Normal file
View File

@@ -0,0 +1,27 @@
function obj = partition(obj)
arguments (Input)
obj (1, 1) {mustBeA(obj, 'miSim')};
end
arguments (Output)
obj (1, 1) {mustBeA(obj, 'miSim')};
end
% Assess sensing performance of each agent at each sample point
% in the domain
agentPerformances = cellfun(@(x) reshape(x.sensorModel.sensorPerformance(x.pos, x.pan, x.tilt, [obj.objective.X(:), obj.objective.Y(:), zeros(size(obj.objective.X(:)))]), size(obj.objective.X)), obj.agents, 'UniformOutput', false);
agentPerformances{end + 1} = obj.sensorPerformanceMinimum * ones(size(agentPerformances{end})); % add additional layer to represent the threshold that has to be cleared for assignment to any partiton
agentPerformances = cat(3, agentPerformances{:});
% Get highest performance value at each point
[~, idx] = max(agentPerformances, [], 3);
% Collect agent indices in the same way
agentInds = cellfun(@(x) x.index * ones(size(obj.objective.X)), obj.agents, 'UniformOutput', false);
agentInds{end + 1} = zeros(size(agentInds{end})); % index for no assignment
agentInds = cat(3, agentInds{:});
% Get highest performing agent's index
[m,n,~] = size(agentInds);
[i,j] = ndgrid(1:m, 1:n);
obj.partitioning = agentInds(sub2ind(size(agentInds), i, j, idx));
end

41
@miSim/plot.m Normal file
View File

@@ -0,0 +1,41 @@
function [obj, f] = plot(obj)
arguments (Input)
obj (1, 1) {mustBeA(obj, 'miSim')};
end
arguments (Output)
obj (1, 1) {mustBeA(obj, 'miSim')};
f (1, 1) {mustBeA(f, 'matlab.ui.Figure')};
end
% Plot domain
[obj.domain, f] = obj.domain.plotWireframe(obj.spatialPlotIndices);
% Plot obstacles
for ii = 1:size(obj.obstacles, 1)
[obj.obstacles{ii}, f] = obj.obstacles{ii}.plotWireframe(obj.spatialPlotIndices, f);
end
% Plot objective gradient
f = obj.domain.objective.plot(obj.objectivePlotIndices, f);
% Plot agents and their collision geometries
for ii = 1:size(obj.agents, 1)
[obj.agents{ii}, f] = obj.agents{ii}.plot(obj.spatialPlotIndices, f);
end
% Plot communication links
[obj, f] = obj.plotConnections(obj.spatialPlotIndices, f);
% Plot abstract network graph
[obj, f] = obj.plotGraph(obj.networkGraphIndex, f);
% Plot domain partitioning
[obj, f] = obj.plotPartitions(obj.partitionGraphIndex, f);
% Enforce plot limits
for ii = 1:size(obj.spatialPlotIndices, 2)
xlim(f.Children(1).Children(obj.spatialPlotIndices(ii)), [obj.domain.minCorner(1), obj.domain.maxCorner(1)]);
ylim(f.Children(1).Children(obj.spatialPlotIndices(ii)), [obj.domain.minCorner(2), obj.domain.maxCorner(2)]);
zlim(f.Children(1).Children(obj.spatialPlotIndices(ii)), [obj.domain.minCorner(3), obj.domain.maxCorner(3)]);
end
end

45
@miSim/plotConnections.m Normal file
View File

@@ -0,0 +1,45 @@
function [obj, f] = plotConnections(obj, ind, f)
arguments (Input)
obj (1, 1) {mustBeA(obj, 'miSim')};
ind (1, :) double = NaN;
f (1, 1) {mustBeA(f, 'matlab.ui.Figure')} = figure;
end
arguments (Output)
obj (1, 1) {mustBeA(obj, 'miSim')};
f (1, 1) {mustBeA(f, 'matlab.ui.Figure')};
end
% Iterate over lower triangle off-diagonal region of the
% adjacency matrix to plot communications links between agents
X = []; Y = []; Z = [];
for ii = 2:size(obj.adjacency, 1)
for jj = 1:(ii - 1)
if obj.adjacency(ii, jj)
X = [X; obj.agents{ii}.pos(1), obj.agents{jj}.pos(1)];
Y = [Y; obj.agents{ii}.pos(2), obj.agents{jj}.pos(2)];
Z = [Z; obj.agents{ii}.pos(3), obj.agents{jj}.pos(3)];
end
end
end
X = X'; Y = Y'; Z = Z';
% Plot the connections
if isnan(ind)
hold(f.CurrentAxes, "on");
o = plot3(f.CurrentAxes, X, Y, Z, 'Color', 'g', 'LineWidth', 2, 'LineStyle', '--');
hold(f.CurrentAxes, "off");
else
hold(f.Children(1).Children(ind(1)), "on");
o = plot3(f.Children(1).Children(ind(1)), X, Y, Z, 'Color', 'g', 'LineWidth', 2, 'LineStyle', '--');
hold(f.Children(1).Children(ind(1)), "off");
end
% Copy to other plots
if size(ind, 2) > 1
for ii = 2:size(ind, 2)
o = [o, copyobj(o(:, 1), f.Children(1).Children(ind(ii)))];
end
end
obj.connectionsPlot = o;
end

31
@miSim/plotGraph.m Normal file
View File

@@ -0,0 +1,31 @@
function [obj, f] = plotGraph(obj, ind, f)
arguments (Input)
obj (1, 1) {mustBeA(obj, 'miSim')};
ind (1, :) double = NaN;
f (1, 1) {mustBeA(f, 'matlab.ui.Figure')} = figure;
end
arguments (Output)
obj (1, 1) {mustBeA(obj, 'miSim')};
f (1, 1) {mustBeA(f, 'matlab.ui.Figure')};
end
% Form graph from adjacency matrix
G = graph(obj.adjacency, 'omitselfloops');
% Plot graph object
if isnan(ind)
hold(f.CurrentAxes, 'on');
o = plot(f.CurrentAxes, G, 'LineStyle', '--', 'EdgeColor', 'g', 'NodeColor', 'k', 'LineWidth', 2);
hold(f.CurrentAxes, 'off');
else
hold(f.Children(1).Children(ind(1)), 'on');
o = plot(f.Children(1).Children(ind(1)), G, 'LineStyle', '--', 'EdgeColor', 'g', 'NodeColor', 'k', 'LineWidth', 2);
hold(f.Children(1).Children(ind(1)), 'off');
if size(ind, 2) > 1
for ii = 2:size(ind, 2)
o = [o; copyobj(o(1), f.Children(1).Children(ind(ii)))];
end
end
end
obj.graphPlot = o;
end

27
@miSim/plotPartitions.m Normal file
View File

@@ -0,0 +1,27 @@
function [obj, f] = plotPartitions(obj, ind, f)
arguments (Input)
obj (1, 1) {mustBeA(obj, 'miSim')};
ind (1, :) double = NaN;
f (1, 1) {mustBeA(f, 'matlab.ui.Figure')} = figure;
end
arguments (Output)
obj (1, 1) {mustBeA(obj, 'miSim')};
f (1, 1) {mustBeA(f, 'matlab.ui.Figure')};
end
if isnan(ind)
hold(f.CurrentAxes, 'on');
o = imagesc(f.CurrentAxes, obj.partitioning);
hold(f.CurrentAxes, 'off');
else
hold(f.Children(1).Children(ind(1)), 'on');
o = imagesc(f.Children(1).Children(ind(1)), obj.partitioning);
hold(f.Children(1).Children(ind(1)), 'on');
if size(ind, 2) > 1
for ii = 2:size(ind, 2)
o = [o, copyobj(o(1), f.Children(1).Children(ind(ii)))];
end
end
end
obj.partitionPlot = o;
end

52
@miSim/run.m Normal file
View File

@@ -0,0 +1,52 @@
function [obj, f] = run(obj, f)
arguments (Input)
obj (1, 1) {mustBeA(obj, 'miSim')};
f (1, 1) {mustBeA(f, 'matlab.ui.Figure')} = figure;
end
arguments (Output)
obj (1, 1) {mustBeA(obj, 'miSim')};
f (1, 1) {mustBeA(f, 'matlab.ui.Figure')};
end
% Create axes if they don't already exist
f = firstPlotSetup(f);
% Set up times to iterate over
times = linspace(0, obj.timestep * obj.maxIter, obj.maxIter+1)';
partitioningTimes = times(obj.partitioningFreq:obj.partitioningFreq:size(times, 1));
% Start video writer
v = obj.setupVideoWriter();
v.open();
for ii = 1:size(times, 1)
% Display current sim time
t = times(ii);
fprintf("Sim Time: %4.2f (%d/%d)\n", t, ii, obj.maxIter)
% Check if it's time for new partitions
updatePartitions = false;
if ismember(t, partitioningTimes)
updatePartitions = true;
obj = obj.partition();
end
% Iterate over agents to simulate their motion
for jj = 1:size(obj.agents, 1)
obj.agents{jj} = obj.agents{jj}.run(obj.objective, obj.domain, obj.partitioning);
end
% Update adjacency matrix
obj = obj.updateAdjacency;
% Update plots
[obj, f] = obj.updatePlots(f, updatePartitions);
% Write frame in to video
I = getframe(f);
v.writeVideo(I);
end
% Close video file
v.close();
end

View File

@@ -1,6 +1,6 @@
function v = setupVideoWriter(timestep) function v = setupVideoWriter(obj)
arguments (Input) arguments (Input)
timestep (1, 1) double; obj (1, 1) {mustBeA(obj, 'miSim')};
end end
arguments (Output) arguments (Output)
v (1, 1) {mustBeA(v, 'VideoWriter')}; v (1, 1) {mustBeA(v, 'VideoWriter')};
@@ -12,6 +12,6 @@ function v = setupVideoWriter(timestep)
v = VideoWriter(fullfile('sandbox', strcat(string(datetime('now'), 'yyyy_MM_dd_HH_mm_ss'), '_miSimHist')), 'Motion JPEG AVI'); v = VideoWriter(fullfile('sandbox', strcat(string(datetime('now'), 'yyyy_MM_dd_HH_mm_ss'), '_miSimHist')), 'Motion JPEG AVI');
end end
v.FrameRate = 1/timestep; v.FrameRate = 1 / obj.timestep;
v.Quality = 90; v.Quality = 90;
end end

32
@miSim/updateAdjacency.m Normal file
View File

@@ -0,0 +1,32 @@
function obj = updateAdjacency(obj)
arguments (Input)
obj (1, 1) {mustBeA(obj, 'miSim')};
end
arguments (Output)
obj (1, 1) {mustBeA(obj, 'miSim')};
end
% Initialize assuming only self-connections
A = logical(eye(size(obj.agents, 1)));
% Check lower triangle off-diagonal connections
for ii = 2:size(A, 1)
for jj = 1:(ii - 1)
if norm(obj.agents{ii}.pos - obj.agents{jj}.pos) <= min([obj.agents{ii}.comRange, obj.agents{jj}.comRange])
% Make sure that obstacles don't obstruct the line
% of sight, breaking the connection
for kk = 1:size(obj.obstacles, 1)
if ~obj.obstacles{kk}.containsLine(obj.agents{ii}.pos, obj.agents{jj}.pos)
A(ii, jj) = true;
end
end
% need extra handling for cases with no obstacles
if isempty(obj.obstacles)
A(ii, jj) = true;
end
end
end
end
obj.adjacency = A | A';
end

43
@miSim/updatePlots.m Normal file
View File

@@ -0,0 +1,43 @@
function [obj, f] = updatePlots(obj, f, updatePartitions)
arguments (Input)
obj (1, 1) {mustBeA(obj, 'miSim')};
f (1, 1) {mustBeA(f, 'matlab.ui.Figure')} = figure;
updatePartitions (1, 1) logical = false;
end
arguments (Output)
obj (1, 1) {mustBeA(obj, 'miSim')};
f (1, 1) {mustBeA(f, 'matlab.ui.Figure')};
end
% Update agent positions, collision geometries
for ii = 1:size(obj.agents, 1)
obj.agents{ii}.updatePlots();
end
% The remaining updates might be possible to do in a clever way
% that moves existing lines instead of clearing and
% re-plotting, which is much better for performance boost
% Update agent connections plot
delete(obj.connectionsPlot);
[obj, f] = obj.plotConnections(obj.spatialPlotIndices, f);
% Update network graph plot
delete(obj.graphPlot);
[obj, f] = obj.plotGraph(obj.networkGraphIndex, f);
% Update partitioning plot
if updatePartitions
delete(obj.partitionPlot);
[obj, f] = obj.plotPartitions(obj.partitionGraphIndex, f);
end
% reset plot limits to fit domain
for ii = 1:size(obj.spatialPlotIndices, 2)
xlim(f.Children(1).Children(obj.spatialPlotIndices(ii)), [obj.domain.minCorner(1), obj.domain.maxCorner(1)]);
ylim(f.Children(1).Children(obj.spatialPlotIndices(ii)), [obj.domain.minCorner(2), obj.domain.maxCorner(2)]);
zlim(f.Children(1).Children(obj.spatialPlotIndices(ii)), [obj.domain.minCorner(3), obj.domain.maxCorner(3)]);
end
drawnow;
end

View File

@@ -0,0 +1,37 @@
function obj = initialize(obj, objectiveFunction, domain, discretizationStep, protectedRange)
arguments (Input)
obj (1,1) {mustBeA(obj, 'sensingObjective')};
objectiveFunction (1, 1) {mustBeA(objectiveFunction, 'function_handle')};
domain (1, 1) {mustBeGeometry};
discretizationStep (1, 1) double = 1;
protectedRange (1, 1) double = 1;
end
arguments (Output)
obj (1,1) {mustBeA(obj, 'sensingObjective')};
end
obj.groundAlt = domain.minCorner(3);
obj.protectedRange = protectedRange;
% Extract footprint limits
xMin = min(domain.footprint(:, 1));
xMax = max(domain.footprint(:, 1));
yMin = min(domain.footprint(:, 2));
yMax = max(domain.footprint(:, 2));
xGrid = unique([xMin:discretizationStep:xMax, xMax]);
yGrid = unique([yMin:discretizationStep:yMax, yMax]);
% Store grid points for plotting later
[obj.X, obj.Y] = meshgrid(xGrid, yGrid);
% Evaluate function over grid points
obj.objectiveFunction = objectiveFunction;
obj.values = reshape(obj.objectiveFunction(obj.X, obj.Y), size(obj.X));
% store ground position
idx = obj.values == max(obj.values, [], "all");
obj.groundPos = [obj.X(idx), obj.Y(idx)];
assert(domain.distance([obj.groundPos, domain.center(3)]) > protectedRange, "Domain is crowding the sensing objective")
end

View File

@@ -0,0 +1,27 @@
function obj = initializeRandomMvnpdf(obj, domain, discretizationStep, protectedRange)
arguments (Input)
obj (1, 1) {mustBeA(obj, 'sensingObjective')};
domain (1, 1) {mustBeGeometry};
discretizationStep (1, 1) double = 1;
protectedRange (1, 1) double = 1;
end
arguments (Output)
obj (1, 1) {mustBeA(obj, 'sensingObjective')};
end
% Set random objective position
mu = domain.minCorner;
while domain.distance(mu) < protectedRange
mu = domain.random();
end
mu = mu(1:2);
% Set random distribution parameters
sig = [2 + rand * 2, 1; 1, 2 + rand * 2];
% Set up random bivariate normal distribution function
objectiveFunction = @(x, y) mvnpdf([x(:), y(:)], mu, sig);
% Regular initialization
obj = obj.initialize(objectiveFunction, domain, discretizationStep, protectedRange);
end

35
@sensingObjective/plot.m Normal file
View File

@@ -0,0 +1,35 @@
function f = plot(obj, ind, f)
arguments (Input)
obj (1,1) {mustBeA(obj, 'sensingObjective')};
ind (1, :) double = NaN;
f (1,1) {mustBeA(f, 'matlab.ui.Figure')} = figure;
end
arguments (Output)
f (1,1) {mustBeA(f, 'matlab.ui.Figure')};
end
% Create axes if they don't already exist
f = firstPlotSetup(f);
% Plot gradient on the "floor" of the domain
if isnan(ind)
hold(f.CurrentAxes, "on");
o = surf(f.CurrentAxes, obj.X, obj.Y, repmat(obj.groundAlt, size(obj.X)), obj.values ./ max(obj.values, [], "all"), 'EdgeColor', 'none');
o.HitTest = 'off';
o.PickableParts = 'none';
hold(f.CurrentAxes, "off");
else
hold(f.Children(1).Children(ind(1)), "on");
o = surf(f.Children(1).Children(ind(1)), obj.X, obj.Y, repmat(obj.groundAlt, size(obj.X)), obj.values ./ max(obj.values, [], "all"), 'EdgeColor', 'none');
o.HitTest = 'off';
o.PickableParts = 'none';
hold(f.Children(1).Children(ind(1)), "off");
end
% Add to other perspectives
if size(ind, 2) > 1
for ii = 2:size(ind, 2)
copyobj(o, f.Children(1).Children(ind(ii)));
end
end
end

View File

@@ -0,0 +1,20 @@
classdef sensingObjective
% Sensing objective definition parent class
properties (SetAccess = private, GetAccess = public)
label = "";
groundAlt = 0;
groundPos = [0, 0];
discretizationStep = 1;
objectiveFunction = @(x, y) 0; % define objective functions over a grid in this manner
X = [];
Y = [];
values = [];
protectedRange = 1; % keep obstacles from crowding objective
end
methods (Access = public)
[obj] = initialize(obj, objectiveFunction, domain, discretizationStep, protectedRange);
[obj] = initializeRandomMvnpdf(obj, domain, protectedRange, discretizationStep, protectedRange);
[f ] = plot(obj, ind, f);
end
end

167
agent.m
View File

@@ -1,167 +0,0 @@
classdef agent
properties (SetAccess = private, GetAccess = public)
% Identifiers
index = NaN;
label = "";
% Sensor
sensorModel;
sensingLength = 0.05; % length parameter used by sensing function
% Guidance
guidanceModel;
% State
lastPos = NaN(1, 3); % position from previous timestep
pos = NaN(1, 3); % current position
vel = NaN(1, 3); % current velocity
pan = NaN; % pan angle
tilt = NaN; % tilt angle
% Collision
collisionGeometry;
% FOV cone
fovGeometry;
% Communication
comRange = NaN;
% Plotting
scatterPoints;
end
methods (Access = public)
function obj = initialize(obj, pos, vel, pan, tilt, collisionGeometry, sensorModel, guidanceModel, comRange, index, label)
arguments (Input)
obj (1, 1) {mustBeA(obj, 'agent')};
pos (1, 3) double;
vel (1, 3) double;
pan (1, 1) double;
tilt (1, 1) double;
collisionGeometry (1, 1) {mustBeGeometry};
sensorModel (1, 1) {mustBeSensor}
guidanceModel (1, 1) {mustBeA(guidanceModel, 'function_handle')};
comRange (1, 1) double = NaN;
index (1, 1) double = NaN;
label (1, 1) string = "";
end
arguments (Output)
obj (1, 1) {mustBeA(obj, 'agent')};
end
obj.pos = pos;
obj.vel = vel;
obj.pan = pan;
obj.tilt = tilt;
obj.collisionGeometry = collisionGeometry;
obj.sensorModel = sensorModel;
obj.guidanceModel = guidanceModel;
obj.comRange = comRange;
obj.index = index;
obj.label = label;
% Initialize FOV cone
obj.fovGeometry = cone;
obj.fovGeometry = obj.fovGeometry.initialize([obj.pos(1:2), 0], obj.sensorModel.r, obj.pos(3), REGION_TYPE.FOV, sprintf("%s FOV", obj.label));
end
function obj = run(obj, sensingObjective, domain, partitioning)
arguments (Input)
obj (1, 1) {mustBeA(obj, 'agent')};
sensingObjective (1, 1) {mustBeA(sensingObjective, 'sensingObjective')};
domain (1, 1) {mustBeGeometry};
partitioning (:, :) double;
end
arguments (Output)
obj (1, 1) {mustBeA(obj, 'agent')};
end
% Do sensing
[sensedValues, sensedPositions] = obj.sensorModel.sense(obj, sensingObjective, domain, partitioning);
% Determine next planned position
nextPos = obj.guidanceModel(sensedValues, sensedPositions, obj.pos);
% Move to next position
% (dynamics not modeled at this time)
obj.lastPos = obj.pos;
obj.pos = nextPos;
% Calculate movement
d = obj.pos - obj.collisionGeometry.center;
% Reinitialize collision geometry in the new position
obj.collisionGeometry = obj.collisionGeometry.initialize([obj.collisionGeometry.minCorner; obj.collisionGeometry.maxCorner] + d, obj.collisionGeometry.tag, obj.collisionGeometry.label);
end
function updatePlots(obj)
arguments (Input)
obj (1, 1) {mustBeA(obj, 'agent')};
end
arguments (Output)
end
% Scatterplot point positions
for ii = 1:size(obj.scatterPoints, 1)
obj.scatterPoints(ii).XData = obj.pos(1);
obj.scatterPoints(ii).YData = obj.pos(2);
obj.scatterPoints(ii).ZData = obj.pos(3);
end
% Find change in agent position since last timestep
deltaPos = obj.pos - obj.lastPos;
% Collision geometry edges
for jj = 1:size(obj.collisionGeometry.lines, 2)
% Update plotting
for ii = 1:size(obj.collisionGeometry.lines(:, jj), 1)
obj.collisionGeometry.lines(ii, jj).XData = obj.collisionGeometry.lines(ii, jj).XData + deltaPos(1);
obj.collisionGeometry.lines(ii, jj).YData = obj.collisionGeometry.lines(ii, jj).YData + deltaPos(2);
obj.collisionGeometry.lines(ii, jj).ZData = obj.collisionGeometry.lines(ii, jj).ZData + deltaPos(3);
end
end
% Update FOV geometry surfaces
for jj = 1:size(obj.fovGeometry.surface, 2)
% Update each plot
obj.fovGeometry.surface(jj).XData = obj.fovGeometry.surface(jj).XData + deltaPos(1);
obj.fovGeometry.surface(jj).YData = obj.fovGeometry.surface(jj).YData + deltaPos(2);
obj.fovGeometry.surface(jj).ZData = obj.fovGeometry.surface(jj).ZData + deltaPos(3);
end
end
function [obj, f] = plot(obj, ind, f)
arguments (Input)
obj (1, 1) {mustBeA(obj, 'agent')};
ind (1, :) double = NaN;
f (1, 1) {mustBeA(f, 'matlab.ui.Figure')} = figure;
end
arguments (Output)
obj (1, 1) {mustBeA(obj, 'agent')};
f (1, 1) {mustBeA(f, 'matlab.ui.Figure')};
end
% Create axes if they don't already exist
f = firstPlotSetup(f);
% Plot points representing the agent position
hold(f.Children(1).Children(end), "on");
o = scatter3(f.Children(1).Children(end), obj.pos(1), obj.pos(2), obj.pos(3), 'filled', 'ko', 'SizeData', 25);
hold(f.Children(1).Children(end), "off");
% Check if this is a tiled layout figure
if strcmp(f.Children(1).Type, 'tiledlayout')
% Add to other perspectives
o = [o; copyobj(o(1), f.Children(1).Children(2))];
o = [o; copyobj(o(1), f.Children(1).Children(3))];
o = [o; copyobj(o(1), f.Children(1).Children(4))];
end
obj.scatterPoints = o;
% Plot collision geometry
[obj.collisionGeometry, f] = obj.collisionGeometry.plotWireframe(ind, f);
% Plot FOV geometry
[obj.fovGeometry, f] = obj.fovGeometry.plot(ind, f);
end
end
end

22
geometries/@cone/cone.m Normal file
View File

@@ -0,0 +1,22 @@
classdef cone
% Conical geometry
properties (SetAccess = private, GetAccess = public)
% Meta
tag = REGION_TYPE.INVALID;
label = "";
% Spatial
center = NaN;
radius = NaN;
height = NaN;
% Plotting
surface;
n = 32;
end
methods (Access = public)
[obj ] = initialize(obj, center, radius, height, tag, label);
[obj, f] = plot(obj, ind, f);
end
end

View File

@@ -0,0 +1,19 @@
function obj = initialize(obj, center, radius, height, tag, label)
arguments (Input)
obj (1, 1) {mustBeA(obj, 'cone')};
center (1, 3) double;
radius (1, 1) double;
height (1, 1) double;
tag (1, 1) REGION_TYPE = REGION_TYPE.INVALID;
label (1, 1) string = "";
end
arguments (Output)
obj (1, 1) {mustBeA(obj, 'cone')};
end
obj.center = center;
obj.radius = radius;
obj.height = height;
obj.tag = tag;
obj.label = label;
end

43
geometries/@cone/plot.m Normal file
View File

@@ -0,0 +1,43 @@
function [obj, f] = plot(obj, ind, f)
arguments (Input)
obj (1, 1) {mustBeA(obj, 'cone')};
ind (1, :) double = NaN;
f (1, 1) {mustBeA(f, 'matlab.ui.Figure')} = figure;
end
arguments (Output)
obj (1, 1) {mustBeA(obj, 'cone')};
f (1, 1) {mustBeA(f, 'matlab.ui.Figure')};
end
% Create axes if they don't already exist
f = firstPlotSetup(f);
% Plot cone
[X, Y, Z] = cylinder([obj.radius, 0], obj.n);
% Scale to match height
Z = Z * obj.height;
% Move to center location
X = X + obj.center(1);
Y = Y + obj.center(2);
Z = Z + obj.center(3);
% Plot
if isnan(ind)
o = surf(f.CurrentAxes, X, Y, Z);
else
hold(f.Children(1).Children(ind(1)), "on");
o = surf(f.Children(1).Children(ind(1)), X, Y, Z, ones([size(Z), 1]) .* reshape(obj.tag.color, 1, 1, 3), 'FaceAlpha', 0.25, 'EdgeColor', 'none');
hold(f.Children(1).Children(ind(1)), "off");
end
% Copy to other requested tiles
if numel(ind) > 1
for ii = 2:size(ind, 2)
o = [o, copyobj(o(:, 1), f.Children(1).Children(ind(ii)))];
end
end
obj.surface = o;
end

View File

@@ -0,0 +1,10 @@
function c = contains(obj, pos)
arguments (Input)
obj (1, 1) {mustBeA(obj, 'rectangularPrism')};
pos (:, 3) double;
end
arguments (Output)
c (:, 1) logical
end
c = all(pos >= repmat(obj.minCorner, size(pos, 1), 1), 2) & all(pos <= repmat(obj.maxCorner, size(pos, 1), 1), 2);
end

View File

@@ -0,0 +1,41 @@
function c = containsLine(obj, pos1, pos2)
arguments (Input)
obj (1, 1) {mustBeA(obj, 'rectangularPrism')};
pos1 (1, 3) double;
pos2 (1, 3) double;
end
arguments (Output)
c (1, 1) logical
end
d = pos2 - pos1;
% edge case where the line is parallel to the geometry
if abs(d) < 1e-12
% check if it happens to start or end inside or outside of
% the geometry
if obj.contains(pos1) || obj.contains(pos2)
c = true;
else
c = false;
end
return;
end
tmin = -inf;
tmax = inf;
% Standard case
for ii = 1:3
t1 = (obj.minCorner(ii) - pos1(ii)) / d(ii);
t2 = (obj.maxCorner(ii) - pos2(ii)) / d(ii);
tmin = max(tmin, min(t1, t2));
tmax = min(tmax, max(t1, t2));
if tmin > tmax
c = false;
return;
end
end
c = (tmax >= 0) && (tmin <= 1);
end

View File

@@ -0,0 +1,32 @@
function d = distance(obj, pos)
arguments (Input)
obj (1, 1) {mustBeA(obj, 'rectangularPrism')};
pos (:, 3) double;
end
arguments (Output)
d (:, 1) double
end
if obj.contains(pos)
% Queried point is inside geometry
% find minimum distance to any face
d = min([pos(1) - obj.minCorner(1), ...
pos(2) - obj.minCorner(2), ...
pos(3) - obj.minCorner(3), ...
obj.maxCorner(1) - pos(1), ...
obj.maxCorner(2) - pos(2), ...
obj.maxCorner(3) - pos(3)]);
else
% Queried point is outside geometry
cPos = NaN(1, 3);
for ii = 1:3
if pos(ii) < obj.minCorner(ii)
cPos(ii) = obj.minCorner(ii);
elseif pos(ii) > obj.maxCorner(ii)
cPos(ii) = obj.maxCorner(ii);
else
cPos(ii) = pos(ii);
end
end
d = norm(cPos - pos);
end
end

View File

@@ -0,0 +1,47 @@
function obj = initialize(obj, bounds, tag, label, objectiveFunction, discretizationStep)
arguments (Input)
obj (1, 1) {mustBeA(obj, 'rectangularPrism')};
bounds (2, 3) double;
tag (1, 1) REGION_TYPE = REGION_TYPE.INVALID;
label (1, 1) string = "";
objectiveFunction (1, 1) function_handle = @(x, y) 1;
discretizationStep (1, 1) double = 1;
end
arguments (Output)
obj (1, 1) {mustBeA(obj, 'rectangularPrism')};
end
obj.tag = tag;
obj.label = label;
% Define geometry bounds by LL corner and UR corner
obj.minCorner = bounds(1, 1:3);
obj.maxCorner = bounds(2, 1:3);
% Compute L, W, H
obj.dimensions = [obj.maxCorner(1) - obj.minCorner(1), obj.maxCorner(2) - obj.minCorner(2), obj.maxCorner(3) - obj.minCorner(3)];
% Compute center
obj.center = obj.minCorner + obj.dimensions ./ 2;
% Compute vertices
obj.vertices = [obj.minCorner;
obj.maxCorner(1), obj.minCorner(2:3);
obj.maxCorner(1:2), obj.minCorner(3);
obj.minCorner(1), obj.maxCorner(2), obj.minCorner(3);
obj.minCorner(1:2), obj.maxCorner(3);
obj.maxCorner(1), obj.minCorner(2), obj.maxCorner(3);
obj.minCorner(1), obj.maxCorner(2:3)
obj.maxCorner;];
% Compute footprint
obj.footprint = [obj.minCorner(1:2); ...
[obj.minCorner(1), obj.maxCorner(2)]; ...
[obj.maxCorner(1), obj.minCorner(2)]; ...
obj.maxCorner(1:2)];
% Instantiate sensingObjective only for DOMAIN-type regions
if tag == REGION_TYPE.DOMAIN
obj.objective = sensingObjective;
end
end

View File

@@ -0,0 +1,44 @@
function [obj] = initializeRandom(obj, tag, label, minDimension, maxDimension, domain)
arguments (Input)
obj (1, 1) {mustBeA(obj, 'rectangularPrism')};
tag (1, 1) REGION_TYPE = REGION_TYPE.INVALID;
label (1, 1) string = "";
minDimension (1, 1) double = 10;
maxDimension (1, 1) double= 20;
domain (1, 1) {mustBeGeometry} = rectangularPrism;
end
arguments (Output)
obj (1, 1) {mustBeA(obj, 'rectangularPrism')};
end
% Produce random bounds based on region type
if tag == REGION_TYPE.DOMAIN
% Domain
L = ceil(minDimension + rand * (maxDimension - minDimension));
bounds = [zeros(1, 3); L * ones(1, 3)];
else
% Obstacle
% Produce a corners that are contained in the domain
ii = 0;
candidateMaxCorner = domain.maxCorner + ones(1, 3);
candidateMinCorner = domain.minCorner - ones(1, 3);
% Continue until the domain contains the obstacle without crowding the objective
while ~domain.contains(candidateMaxCorner) || all(domain.objective.groundPos + domain.objective.protectedRange >= candidateMinCorner(1:2), 2) && all(domain.objective.groundPos - domain.objective.protectedRange <= candidateMaxCorner(1:2), 2)
if ii == 0 || ii > 10
candidateMinCorner = domain.random();
candidateMinCorner(3) = 0; % bind to floor
ii = 1;
end
candidateMaxCorner = candidateMinCorner + minDimension + rand(1, 3) * (maxDimension - minDimension);
ii = ii + 1;
end
bounds = [candidateMinCorner; candidateMaxCorner;];
end
% Regular initialization
obj = obj.initialize(bounds, tag, label);
end

View File

@@ -0,0 +1,37 @@
function [obj, f] = plotWireframe(obj, ind, f)
arguments (Input)
obj (1, 1) {mustBeA(obj, 'rectangularPrism')};
ind (1, :) double = NaN;
f (1, 1) {mustBeA(f, 'matlab.ui.Figure')} = figure;
end
arguments (Output)
obj (1, 1) {mustBeA(obj, 'rectangularPrism')};
f (1, 1) {mustBeA(f, 'matlab.ui.Figure')};
end
% Create axes if they don't already exist
f = firstPlotSetup(f);
% Create plotting inputs from vertices and edges
X = [obj.vertices(obj.edges(:,1),1), obj.vertices(obj.edges(:,2),1)]';
Y = [obj.vertices(obj.edges(:,1),2), obj.vertices(obj.edges(:,2),2)]';
Z = [obj.vertices(obj.edges(:,1),3), obj.vertices(obj.edges(:,2),3)]';
% Plot the boundaries of the geometry into 3D view
if isnan(ind)
o = plot3(f.CurrentAxes, X, Y, Z, '-', 'Color', obj.tag.color, 'LineWidth', 2);
else
hold(f.Children(1).Children(ind(1)), "on");
o = plot3(f.Children(1).Children(ind(1)), X, Y, Z, '-', 'Color', obj.tag.color, 'LineWidth', 2);
hold(f.Children(1).Children(ind(1)), "off");
end
% Copy to other requested tiles
if numel(ind) > 1
for ii = 2:size(ind, 2)
o = [o, copyobj(o(:, 1), f.Children(1).Children(ind(ii)))];
end
end
obj.lines = o;
end

View File

@@ -0,0 +1,9 @@
function r = random(obj)
arguments (Input)
obj (1, 1) {mustBeA(obj, 'rectangularPrism')};
end
arguments (Output)
r (1, 3) double
end
r = (obj.vertices(1, 1:3) + rand(1, 3) .* obj.vertices(8, 1:3) - obj.vertices(1, 1:3))';
end

View File

@@ -0,0 +1,38 @@
classdef rectangularPrism
% Rectangular prism geometry
properties (SetAccess = private, GetAccess = public)
% Meta
tag = REGION_TYPE.INVALID;
label = "";
% Spatial
minCorner = NaN(1, 3);
maxCorner = NaN(1, 3);
dimensions = NaN(1, 3);
center = NaN;
footprint = NaN(4, 2);
% Graph
vertices = NaN(8, 3);
edges = [1 2; 2 3; 3 4; 4 1; % bottom square
5 6; 6 8; 8 7; 7 5; % top square
1 5; 2 6; 3 8; 4 7]; % vertical edges
% Plotting
lines;
end
properties (SetAccess = public, GetAccess = public)
% Sensing objective (for DOMAIN region type only)
objective;
end
methods (Access = public)
[obj ] = initialize(obj, bounds, tag, label, objectiveFunction, discretizationStep);
[obj ] = initializeRandom(obj, tag, label, minDimension, maxDimension, domain);
[r ] = random(obj);
[c ] = contains(obj, pos);
[d ] = distance(obj, pos);
[c ] = containsLine(obj, pos1, pos2);
[obj, f] = plotWireframe(obj, ind, f);
end
end

View File

@@ -1,82 +0,0 @@
classdef cone
% Conical geometry
properties (SetAccess = private, GetAccess = public)
% Meta
tag = REGION_TYPE.INVALID;
label = "";
% Spatial
center = NaN;
radius = NaN;
height = NaN;
% Plotting
surface;
n = 32;
end
methods
function obj = initialize(obj, center, radius, height, tag, label)
arguments (Input)
obj (1, 1) {mustBeA(obj, 'cone')};
center (1, 3) double;
radius (1, 1) double;
height (1, 1) double;
tag (1, 1) REGION_TYPE = REGION_TYPE.INVALID;
label (1, 1) string = "";
end
arguments (Output)
obj (1, 1) {mustBeA(obj, 'cone')};
end
obj.center = center;
obj.radius = radius;
obj.height = height;
obj.tag = tag;
obj.label = label;
end
function [obj, f] = plot(obj, ind, f)
arguments (Input)
obj (1, 1) {mustBeA(obj, 'cone')};
ind (1, :) double = NaN;
f (1, 1) {mustBeA(f, 'matlab.ui.Figure')} = figure;
end
arguments (Output)
obj (1, 1) {mustBeA(obj, 'cone')};
f (1, 1) {mustBeA(f, 'matlab.ui.Figure')};
end
% Create axes if they don't already exist
f = firstPlotSetup(f);
% Plot cone
[X, Y, Z] = cylinder([obj.radius, 0], obj.n);
% Scale to match height
Z = Z * obj.height;
% Move to center location
X = X + obj.center(1);
Y = Y + obj.center(2);
Z = Z + obj.center(3);
% Plot
if isnan(ind)
o = surf(f.CurrentAxes, X, Y, Z);
else
hold(f.Children(1).Children(ind(1)), "on");
o = surf(f.Children(1).Children(ind(1)), X, Y, Z, ones([size(Z), 1]) .* reshape(obj.tag.color, 1, 1, 3), 'FaceAlpha', 0.25, 'EdgeColor', 'none');
hold(f.Children(1).Children(ind(1)), "off");
end
% Copy to other requested tiles
if numel(ind) > 1
for ii = 2:size(ind, 2)
o = [o, copyobj(o(:, 1), f.Children(1).Children(ind(ii)))];
end
end
obj.surface = o;
end
end
end

View File

@@ -1,204 +0,0 @@
classdef rectangularPrism
% Rectangular prism geometry
properties (SetAccess = private, GetAccess = public)
% Meta
tag = REGION_TYPE.INVALID;
label = "";
% Spatial
minCorner = NaN(1, 3);
maxCorner = NaN(1, 3);
dimensions = NaN(1, 3);
center = NaN;
footprint = NaN(4, 2);
% Graph
vertices = NaN(8, 3);
edges = [1 2; 2 3; 3 4; 4 1; % bottom square
5 6; 6 8; 8 7; 7 5; % top square
1 5; 2 6; 3 8; 4 7]; % vertical edges
% Plotting
lines;
end
methods (Access = public)
function obj = initialize(obj, bounds, tag, label)
arguments (Input)
obj (1, 1) {mustBeA(obj, 'rectangularPrism')};
bounds (2, 3) double;
tag (1, 1) REGION_TYPE = REGION_TYPE.INVALID;
label (1, 1) string = "";
end
arguments (Output)
obj (1, 1) {mustBeA(obj, 'rectangularPrism')};
end
obj.tag = tag;
obj.label = label;
%% Define geometry bounds by LL corner and UR corner
obj.minCorner = bounds(1, 1:3);
obj.maxCorner = bounds(2, 1:3);
% Compute L, W, H
obj.dimensions = [obj.maxCorner(1) - obj.minCorner(1), obj.maxCorner(2) - obj.minCorner(2), obj.maxCorner(3) - obj.minCorner(3)];
% Compute center
obj.center = obj.minCorner + obj.dimensions ./ 2;
% Compute vertices
obj.vertices = [obj.minCorner;
obj.maxCorner(1), obj.minCorner(2:3);
obj.maxCorner(1:2), obj.minCorner(3);
obj.minCorner(1), obj.maxCorner(2), obj.minCorner(3);
obj.minCorner(1:2), obj.maxCorner(3);
obj.maxCorner(1), obj.minCorner(2), obj.maxCorner(3);
obj.minCorner(1), obj.maxCorner(2:3)
obj.maxCorner;];
% Compute footprint
obj.footprint = [obj.minCorner(1:2); ...
[obj.minCorner(1), obj.maxCorner(2)]; ...
[obj.maxCorner(1), obj.minCorner(2)]; ...
obj.maxCorner(1:2)];
end
function r = random(obj)
arguments (Input)
obj (1, 1) {mustBeA(obj, 'rectangularPrism')};
end
arguments (Output)
r (1, 3) double
end
r = (obj.vertices(1, 1:3) + rand(1, 3) .* obj.vertices(8, 1:3) - obj.vertices(1, 1:3))';
end
function d = distance(obj, pos)
arguments (Input)
obj (1, 1) {mustBeA(obj, 'rectangularPrism')};
pos (:, 3) double;
end
arguments (Output)
d (:, 1) double
end
assert(~obj.contains(pos), "Cannot determine distance for a point inside of the geometry");
cPos = NaN(1, 3);
for ii = 1:3
if pos(ii) < obj.minCorner(ii)
cPos(ii) = obj.minCorner(ii);
elseif pos(ii) > obj.maxCorner(ii)
cPos(ii) = obj.maxCorner(ii);
else
cPos(ii) = pos(ii);
end
end
d = norm(cPos - pos);
end
function d = interiorDistance(obj, pos)
arguments (Input)
obj (1, 1) {mustBeA(obj, 'rectangularPrism')};
pos (:, 3) double;
end
arguments (Output)
d (:, 1) double
end
assert(obj.contains(pos), "Cannot determine interior distance for a point outside of the geometry");
% find minimum distance to any face
d = min([pos(1) - obj.minCorner(1), ...
pos(2) - obj.minCorner(2), ...
pos(3) - obj.minCorner(3), ...
obj.maxCorner(1) - pos(1), ...
obj.maxCorner(2) - pos(2), ...
obj.maxCorner(3) - pos(3)]);
end
function c = contains(obj, pos)
arguments (Input)
obj (1, 1) {mustBeA(obj, 'rectangularPrism')};
pos (:, 3) double;
end
arguments (Output)
c (:, 1) logical
end
c = all(pos >= repmat(obj.minCorner, size(pos, 1), 1), 2) & all(pos <= repmat(obj.maxCorner, size(pos, 1), 1), 2);
end
function c = containsLine(obj, pos1, pos2)
arguments (Input)
obj (1, 1) {mustBeA(obj, 'rectangularPrism')};
pos1 (1, 3) double;
pos2 (1, 3) double;
end
arguments (Output)
c (1, 1) logical
end
d = pos2 - pos1;
% edge case where the line is parallel to the geometry
if abs(d) < 1e-12
% check if it happens to start or end inside or outside of
% the geometry
if obj.contains(pos1) || obj.contains(pos2)
c = true;
else
c = false;
end
return;
end
tmin = -inf;
tmax = inf;
% Standard case
for ii = 1:3
t1 = (obj.minCorner(ii) - pos1(ii)) / d(ii);
t2 = (obj.maxCorner(ii) - pos2(ii)) / d(ii);
tmin = max(tmin, min(t1, t2));
tmax = min(tmax, max(t1, t2));
if tmin > tmax
c = false;
return;
end
end
c = (tmax >= 0) && (tmin <= 1);
end
function [obj, f] = plotWireframe(obj, ind, f)
arguments (Input)
obj (1, 1) {mustBeA(obj, 'rectangularPrism')};
ind (1, :) double = NaN;
f (1, 1) {mustBeA(f, 'matlab.ui.Figure')} = figure;
end
arguments (Output)
obj (1, 1) {mustBeA(obj, 'rectangularPrism')};
f (1, 1) {mustBeA(f, 'matlab.ui.Figure')};
end
% Create axes if they don't already exist
f = firstPlotSetup(f);
% Create plotting inputs from vertices and edges
X = [obj.vertices(obj.edges(:,1),1), obj.vertices(obj.edges(:,2),1)]';
Y = [obj.vertices(obj.edges(:,1),2), obj.vertices(obj.edges(:,2),2)]';
Z = [obj.vertices(obj.edges(:,1),3), obj.vertices(obj.edges(:,2),3)]';
% Plot the boundaries of the geometry into 3D view
if isnan(ind)
o = plot3(f.CurrentAxes, X, Y, Z, '-', 'Color', obj.tag.color, 'LineWidth', 2);
else
hold(f.Children(1).Children(ind(1)), "on");
o = plot3(f.Children(1).Children(ind(1)), X, Y, Z, '-', 'Color', obj.tag.color, 'LineWidth', 2);
hold(f.Children(1).Children(ind(1)), "off");
end
% Copy to other requested tiles
if numel(ind) > 1
for ii = 2:size(ind, 2)
o = [o, copyobj(o(:, 1), f.Children(1).Children(ind(ii)))];
end
end
obj.lines = o;
end
end
end

379
miSim.m
View File

@@ -1,379 +0,0 @@
classdef miSim
% multiagent interconnection simulation
% Simulation parameters
properties (SetAccess = private, GetAccess = public)
timestep = NaN; % delta time interval for simulation iterations
partitioningFreq = NaN; % number of simulation timesteps at which the partitioning routine is re-run
maxIter = NaN; % maximum number of simulation iterations
domain = rectangularPrism;
objective = sensingObjective;
obstacles = cell(0, 1); % geometries that define obstacles within the domain
agents = cell(0, 1); % agents that move within the domain
adjacency = NaN; % Adjacency matrix representing communications network graph
partitioning = NaN;
end
properties (Access = private)
% Plot objects
connectionsPlot; % objects for lines connecting agents in spatial plots
graphPlot; % objects for abstract network graph plot
partitionPlot; % objects for partition plot
% Indicies for various plot types in the main tiled layout figure
spatialPlotIndices = [6, 4, 3, 2];
objectivePlotIndices = [6, 4];
networkGraphIndex = 5;
partitionGraphIndex = 1;
end
methods (Access = public)
function [obj, f] = initialize(obj, domain, objective, agents, timestep, partitoningFreq, maxIter, obstacles)
arguments (Input)
obj (1, 1) {mustBeA(obj, 'miSim')};
domain (1, 1) {mustBeGeometry};
objective (1, 1) {mustBeA(objective, 'sensingObjective')};
agents (:, 1) cell;
timestep (:, 1) double = 0.05;
partitoningFreq (:, 1) double = 0.25
maxIter (:, 1) double = 1000;
obstacles (:, 1) cell {mustBeGeometry} = cell(0, 1);
end
arguments (Output)
obj (1, 1) {mustBeA(obj, 'miSim')};
f (1, 1) {mustBeA(f, 'matlab.ui.Figure')};
end
% Define simulation time parameters
obj.timestep = timestep;
obj.maxIter = maxIter;
% Define domain
obj.domain = domain;
obj.partitioningFreq = partitoningFreq;
% Add geometries representing obstacles within the domain
obj.obstacles = obstacles;
% Define objective
obj.objective = objective;
% Define agents
obj.agents = agents;
% Compute adjacency matrix
obj = obj.updateAdjacency();
% Create initial partitioning
obj = obj.partition();
% Set up initial plot
% Set up axes arrangement
% Plot domain
[obj.domain, f] = obj.domain.plotWireframe(obj.spatialPlotIndices);
% Plot obstacles
for ii = 1:size(obj.obstacles, 1)
[obj.obstacles{ii}, f] = obj.obstacles{ii}.plotWireframe(obj.spatialPlotIndices, f);
end
% Plot objective gradient
f = obj.objective.plot(obj.objectivePlotIndices, f);
% Plot agents and their collision geometries
for ii = 1:size(obj.agents, 1)
[obj.agents{ii}, f] = obj.agents{ii}.plot(obj.spatialPlotIndices, f);
end
% Plot communication links
[obj, f] = obj.plotConnections(obj.spatialPlotIndices, f);
% Plot abstract network graph
[obj, f] = obj.plotGraph(obj.networkGraphIndex, f);
% Plot domain partitioning
[obj, f] = obj.plotPartitions(obj.partitionGraphIndex, f);
end
function [obj, f] = run(obj, f)
arguments (Input)
obj (1, 1) {mustBeA(obj, 'miSim')};
f (1, 1) {mustBeA(f, 'matlab.ui.Figure')} = figure;
end
arguments (Output)
obj (1, 1) {mustBeA(obj, 'miSim')};
f (1, 1) {mustBeA(f, 'matlab.ui.Figure')};
end
% Create axes if they don't already exist
f = firstPlotSetup(f);
% Set up times to iterate over
times = linspace(0, obj.timestep * obj.maxIter, obj.maxIter+1)';
partitioningTimes = times(obj.partitioningFreq:obj.partitioningFreq:size(times, 1));
% Start video writer
v = setupVideoWriter(obj.timestep);
v.open();
for ii = 1:size(times, 1)
% Display current sim time
t = times(ii);
fprintf("Sim Time: %4.2f (%d/%d)\n", t, ii, obj.maxIter)
% Check if it's time for new partitions
updatePartitions = false;
if ismember(t, partitioningTimes)
updatePartitions = true;
obj = obj.partition();
end
% Iterate over agents to simulate their motion
for jj = 1:size(obj.agents, 1)
obj.agents{jj} = obj.agents{jj}.run(obj.objective, obj.domain, obj.partitioning);
end
% Update adjacency matrix
obj = obj.updateAdjacency;
% Update plots
[obj, f] = obj.updatePlots(f, updatePartitions);
% Write frame in to video
I = getframe(f);
v.writeVideo(I);
end
% Close video file
v.close();
end
function obj = partition(obj)
arguments (Input)
obj (1, 1) {mustBeA(obj, 'miSim')};
end
arguments (Output)
obj (1, 1) {mustBeA(obj, 'miSim')};
end
% Assess sensing performance of each agent at each sample point
% in the domain
agentPerformances = cellfun(@(x) reshape(x.sensorModel.sensorPerformance(x.pos, x.pan, x.tilt, [obj.objective.X(:), obj.objective.Y(:), zeros(size(obj.objective.X(:)))]), size(obj.objective.X)), obj.agents, 'UniformOutput', false);
agentPerformances = cat(3, agentPerformances{:});
% Get highest performance value at each point
[~, idx] = max(agentPerformances, [], 3);
% Collect agent indices in the same way
agentInds = cellfun(@(x) x.index * ones(size(obj.objective.X)), obj.agents, 'UniformOutput', false);
agentInds = cat(3, agentInds{:});
% Get highest performing agent's index
[m,n,~] = size(agentInds);
[i,j] = ndgrid(1:m, 1:n);
obj.partitioning = agentInds(sub2ind(size(agentInds), i, j, idx));
end
function [obj, f] = updatePlots(obj, f, updatePartitions)
arguments (Input)
obj (1, 1) {mustBeA(obj, 'miSim')};
f (1, 1) {mustBeA(f, 'matlab.ui.Figure')} = figure;
updatePartitions (1, 1) logical = false;
end
arguments (Output)
obj (1, 1) {mustBeA(obj, 'miSim')};
f (1, 1) {mustBeA(f, 'matlab.ui.Figure')};
end
% Update agent positions, collision geometries
for ii = 1:size(obj.agents, 1)
obj.agents{ii}.updatePlots();
end
% The remaining updates might be possible to do in a clever way
% that moves existing lines instead of clearing and
% re-plotting, which is much better for performance boost
% Update agent connections plot
delete(obj.connectionsPlot);
[obj, f] = obj.plotConnections(obj.spatialPlotIndices, f);
% Update network graph plot
delete(obj.graphPlot);
[obj, f] = obj.plotGraph(obj.networkGraphIndex, f);
% Update partitioning plot
if updatePartitions
delete(obj.partitionPlot);
[obj, f] = obj.plotPartitions(obj.partitionGraphIndex, f);
end
% reset plot limits to fit domain
for ii = 1:size(obj.spatialPlotIndices, 2)
xlim(f.Children(1).Children(obj.spatialPlotIndices(ii)), [obj.domain.minCorner(1), obj.domain.maxCorner(1)]);
ylim(f.Children(1).Children(obj.spatialPlotIndices(ii)), [obj.domain.minCorner(2), obj.domain.maxCorner(2)]);
zlim(f.Children(1).Children(obj.spatialPlotIndices(ii)), [obj.domain.minCorner(3), obj.domain.maxCorner(3)]);
end
drawnow;
end
function obj = updateAdjacency(obj)
arguments (Input)
obj (1, 1) {mustBeA(obj, 'miSim')};
end
arguments (Output)
obj (1, 1) {mustBeA(obj, 'miSim')};
end
% Initialize assuming only self-connections
A = logical(eye(size(obj.agents, 1)));
% Check lower triangle off-diagonal connections
for ii = 2:size(A, 1)
for jj = 1:(ii - 1)
if norm(obj.agents{ii}.pos - obj.agents{jj}.pos) <= min([obj.agents{ii}.comRange, obj.agents{jj}.comRange])
% Make sure that obstacles don't obstruct the line
% of sight, breaking the connection
for kk = 1:size(obj.obstacles, 1)
if ~obj.obstacles{kk}.containsLine(obj.agents{ii}.pos, obj.agents{jj}.pos)
A(ii, jj) = true;
end
end
% need extra handling for cases with no obstacles
if isempty(obj.obstacles)
A(ii, jj) = true;
end
end
end
end
obj.adjacency = A | A';
end
function [obj, f] = plotConnections(obj, ind, f)
arguments (Input)
obj (1, 1) {mustBeA(obj, 'miSim')};
ind (1, :) double = NaN;
f (1, 1) {mustBeA(f, 'matlab.ui.Figure')} = figure;
end
arguments (Output)
obj (1, 1) {mustBeA(obj, 'miSim')};
f (1, 1) {mustBeA(f, 'matlab.ui.Figure')};
end
% Iterate over lower triangle off-diagonal region of the
% adjacency matrix to plot communications links between agents
X = []; Y = []; Z = [];
for ii = 2:size(obj.adjacency, 1)
for jj = 1:(ii - 1)
if obj.adjacency(ii, jj)
X = [X; obj.agents{ii}.pos(1), obj.agents{jj}.pos(1)];
Y = [Y; obj.agents{ii}.pos(2), obj.agents{jj}.pos(2)];
Z = [Z; obj.agents{ii}.pos(3), obj.agents{jj}.pos(3)];
end
end
end
X = X'; Y = Y'; Z = Z';
% Plot the connections
if isnan(ind)
hold(f.CurrentAxes, "on");
o = plot3(f.CurrentAxes, X, Y, Z, 'Color', 'g', 'LineWidth', 2, 'LineStyle', '--');
hold(f.CurrentAxes, "off");
else
hold(f.Children(1).Children(ind(1)), "on");
o = plot3(f.Children(1).Children(ind(1)), X, Y, Z, 'Color', 'g', 'LineWidth', 2, 'LineStyle', '--');
hold(f.Children(1).Children(ind(1)), "off");
end
% Copy to other plots
if size(ind, 2) > 1
for ii = 2:size(ind, 2)
o = [o, copyobj(o(:, 1), f.Children(1).Children(ind(ii)))];
end
end
obj.connectionsPlot = o;
end
function [obj, f] = plotPartitions(obj, ind, f)
arguments (Input)
obj (1, 1) {mustBeA(obj, 'miSim')};
ind (1, :) double = NaN;
f (1, 1) {mustBeA(f, 'matlab.ui.Figure')} = figure;
end
arguments (Output)
obj (1, 1) {mustBeA(obj, 'miSim')};
f (1, 1) {mustBeA(f, 'matlab.ui.Figure')};
end
if isnan(ind)
hold(f.CurrentAxes, 'on');
o = imagesc(f.CurrentAxes, obj.partitioning);
hold(f.CurrentAxes, 'off');
else
hold(f.Children(1).Children(ind(1)), 'on');
o = imagesc(f.Children(1).Children(ind(1)), obj.partitioning);
hold(f.Children(1).Children(ind(1)), 'on');
if size(ind, 2) > 1
for ii = 2:size(ind, 2)
o = [o, copyobj(o(1), f.Children(1).Children(ind(ii)))];
end
end
end
obj.partitionPlot = o;
end
function [obj, f] = plotGraph(obj, ind, f)
arguments (Input)
obj (1, 1) {mustBeA(obj, 'miSim')};
ind (1, :) double = NaN;
f (1, 1) {mustBeA(f, 'matlab.ui.Figure')} = figure;
end
arguments (Output)
obj (1, 1) {mustBeA(obj, 'miSim')};
f (1, 1) {mustBeA(f, 'matlab.ui.Figure')};
end
% Form graph from adjacency matrix
G = graph(obj.adjacency, 'omitselfloops');
% Plot graph object
if isnan(ind)
hold(f.CurrentAxes, 'on');
o = plot(f.CurrentAxes, G, 'LineStyle', '--', 'EdgeColor', 'g', 'NodeColor', 'k', 'LineWidth', 2);
hold(f.CurrentAxes, 'off');
else
hold(f.Children(1).Children(ind(1)), 'on');
o = plot(f.Children(1).Children(ind(1)), G, 'LineStyle', '--', 'EdgeColor', 'g', 'NodeColor', 'k', 'LineWidth', 2);
hold(f.Children(1).Children(ind(1)), 'off');
if size(ind, 2) > 1
for ii = 2:size(ind, 2)
o = [o; copyobj(o(1), f.Children(1).Children(ind(ii)))];
end
end
end
obj.graphPlot = o;
end
end
methods (Access = private)
function validateInitialization(obj)
% Assert obstacles do not intersect with the domain
% Assert obstacles do not intersect with each other
% Assert the objective has only one maxima within the domain
% Assert the objective's sole maximum is not inaccessible due
% to the placement of an obstacle
end
function validateLoop(obj)
% Assert that agents are safely inside the domain
% Assert that agents are not in proximity to obstacles
% Assert that agents are not in proximity to each other
% Assert that agents form a connected graph
end
end
end

View File

@@ -0,0 +1,2 @@
<?xml version="1.0" encoding="UTF-8"?>
<Info location="sensorPerformance.m" type="File"/>

View File

@@ -0,0 +1,2 @@
<?xml version="1.0" encoding="UTF-8"?>
<Info location="initialize.m" type="File"/>

View File

@@ -0,0 +1,2 @@
<?xml version="1.0" encoding="UTF-8"?>
<Info location="sense.m" type="File"/>

View File

@@ -0,0 +1,2 @@
<?xml version="1.0" encoding="UTF-8"?>
<Info location="1" type="DIR_SIGNIFIER"/>

View File

@@ -0,0 +1,2 @@
<?xml version="1.0" encoding="UTF-8"?>
<Info Ref="util/validators" Type="Relative"/>

View File

@@ -0,0 +1,2 @@
<?xml version="1.0" encoding="UTF-8"?>
<Info location="1bb12f2e-385b-41a4-83e8-f9a9326d95ee" type="Reference"/>

View File

@@ -0,0 +1,2 @@
<?xml version="1.0" encoding="UTF-8"?>
<Info Ref="test" Type="Relative"/>

View File

@@ -0,0 +1,2 @@
<?xml version="1.0" encoding="UTF-8"?>
<Info location="89c37511-fb1f-420e-919a-c5b38c02b501" type="Reference"/>

View File

@@ -0,0 +1,2 @@
<?xml version="1.0" encoding="UTF-8"?>
<Info Ref="util/validators/arguments" Type="Relative"/>

View File

@@ -0,0 +1,2 @@
<?xml version="1.0" encoding="UTF-8"?>
<Info location="30f974f9-50a8-405c-9a83-e7fd84333f0e" type="Reference"/>

View File

@@ -1,2 +0,0 @@
<?xml version="1.0" encoding="UTF-8"?>
<Info Ref="validators" Type="Relative"/>

View File

@@ -1,2 +0,0 @@
<?xml version="1.0" encoding="UTF-8"?>
<Info location="5f96e524-3aac-4fa1-95df-67fd6ce02ff3" type="Reference"/>

View File

@@ -1,2 +0,0 @@
<?xml version="1.0" encoding="UTF-8"?>
<Info Ref="validators/arguments" Type="Relative"/>

View File

@@ -1,2 +0,0 @@
<?xml version="1.0" encoding="UTF-8"?>
<Info location="b7c7eec5-a318-4c17-adb2-b13a21bf0609" type="Reference"/>

View File

@@ -0,0 +1,2 @@
<?xml version="1.0" encoding="UTF-8"?>
<Info Ref="util" Type="Relative"/>

View File

@@ -0,0 +1,2 @@
<?xml version="1.0" encoding="UTF-8"?>
<Info location="2e3f60de-3b82-4ad5-af81-57781353dcbf" type="Reference"/>

View File

@@ -0,0 +1,2 @@
<?xml version="1.0" encoding="UTF-8"?>
<Info/>

View File

@@ -0,0 +1,2 @@
<?xml version="1.0" encoding="UTF-8"?>
<Info location="1" type="DIR_SIGNIFIER"/>

View File

@@ -0,0 +1,2 @@
<?xml version="1.0" encoding="UTF-8"?>
<Info location="plot.m" type="File"/>

View File

@@ -0,0 +1,2 @@
<?xml version="1.0" encoding="UTF-8"?>
<Info location="initialize.m" type="File"/>

View File

@@ -0,0 +1,2 @@
<?xml version="1.0" encoding="UTF-8"?>
<Info/>

View File

@@ -0,0 +1,2 @@
<?xml version="1.0" encoding="UTF-8"?>
<Info location="1" type="DIR_SIGNIFIER"/>

View File

@@ -1,2 +0,0 @@
<?xml version="1.0" encoding="UTF-8"?>
<Info location="mustBeAgents.m" type="File"/>

View File

@@ -1,2 +0,0 @@
<?xml version="1.0" encoding="UTF-8"?>
<Info location="mustBeDcm.m" type="File"/>

View File

@@ -0,0 +1,2 @@
<?xml version="1.0" encoding="UTF-8"?>
<Info location="sense.m" type="File"/>

View File

@@ -0,0 +1,2 @@
<?xml version="1.0" encoding="UTF-8"?>
<Info/>

View File

@@ -0,0 +1,2 @@
<?xml version="1.0" encoding="UTF-8"?>
<Info location="1" type="DIR_SIGNIFIER"/>

View File

@@ -0,0 +1,2 @@
<?xml version="1.0" encoding="UTF-8"?>
<Info location="initialize.m" type="File"/>

View File

@@ -0,0 +1,2 @@
<?xml version="1.0" encoding="UTF-8"?>
<Info location="sensorPerformance.m" type="File"/>

View File

@@ -0,0 +1,2 @@
<?xml version="1.0" encoding="UTF-8"?>
<Info/>

View File

@@ -0,0 +1,2 @@
<?xml version="1.0" encoding="UTF-8"?>
<Info location="@cone" type="File"/>

View File

@@ -0,0 +1,2 @@
<?xml version="1.0" encoding="UTF-8"?>
<Info/>

View File

@@ -0,0 +1,2 @@
<?xml version="1.0" encoding="UTF-8"?>
<Info location="@rectangularPrism" type="File"/>

View File

@@ -0,0 +1,2 @@
<?xml version="1.0" encoding="UTF-8"?>
<Info location="updateAdjacency.m" type="File"/>

View File

@@ -0,0 +1,2 @@
<?xml version="1.0" encoding="UTF-8"?>
<Info/>

Some files were not shown because too many files have changed in this diff Show More