Initial commit

Author: Hespe

.gitignore

@@ -0,0 +1,31 @@
+# Windows default autosave extension
+*.asv
+
+# OSX / *nix default autosave extension
+*.m~
+
+# Compiled MEX binaries (all platforms)
+*.mex*
+
+# Packaged app and toolbox files
+*.mlappinstall
+*.mltbx
+
+# Generated helpsearch folders
+helpsearch*/
+
+# Simulink code generation folders
+slprj/
+sccprj/
+
+# Matlab code generation folders
+codegen/
+
+# Simulink autosave extension
+*.autosave
+
+# Simulink cache files
+*.slxc
+
+# Octave session info
+octave-workspace

conservatism.m

@@ -0,0 +1,75 @@
+clear
+
+addpath('graphs', 'performance')
+
+%% Set up the problem
+kappa = 0.1;
+sysD  = ss(1, 1, 1, 0);
+sysC  = ss(-kappa, 0, 0, 0);
+
+% Grid the probability axis. Higher density where larger changes are
+% expected
+p = [0, 1e-6, 1e-5, 1e-4, 1e-3:1e-3:9e-3, 0.01:0.01:0.29, 0.3:0.025:1];
+
+% Calculate the norms for two examples, large and small networks
+h = [4, 50];
+N = h.*(h+1)/2;
+
+%%
+[H, P] = meshgrid(h, p);
+sz = size(H);
+
+% Pre-allocate storage
+norm_mean  = NaN(length(p),length(h));
+norm_enum  = NaN(length(p),length(h));
+norm_decom = NaN(length(p),length(h));
+
+tic
+parfor i = 1:length(h)*length(p)
+    L = full(laplace_matrix(triangle_graph(H(i))));
+    
+    norm_mean(i)  = performance_mean(sysD, sysC, L, P(i));
+    norm_decom(i) = performance_decomposed(sysD, sysC, L, P(i));
+
+    % The enumerated norm is only calculated for the small network, since
+    % its scalling exponentially and thus infeasible to calculate for the
+    % large network.
+    [~, j] = ind2sub(sz, i);
+    if j == 1
+        norm_enum(i) = performance_enumerated(sysD, sysC, L, P(i), true);
+    end
+end
+toc
+
+%% Plot data like in the paper
+figure(1)
+plot(p, norm_mean(:,1), 'k-.', p, norm_enum(:,1), 'k-', p, norm_decom(:,1), 'k--')
+ylim([0, 44])
+title('Small Network Performance')
+xlabel('Probability p')
+ylabel('H_2-norm')
+legend('Mean', 'Enumerated', 'Decomposed')
+
+figure(2)
+semilogy(p, norm_mean(:,2), 'k-.', p, norm_decom(:,2), 'k--')
+ylim([3e1, 2e4])
+title('Large Network Performance')
+xlabel('Probability p')
+ylabel('H_2-norm')
+legend('Mean', 'Decomposed')
+
+%% CSV Export
+p = p';
+small_mean  = norm_mean(:,1);
+small_enum  = norm_enum(:,1);
+small_decom = norm_decom(:,1);
+small_table = table(p, small_mean, small_enum, small_decom);
+
+large_mean  = norm_mean(:,2);
+large_decom = norm_decom(:,2);
+large_table = table(p, large_mean, large_decom);
+
+name = sprintf('conservatism_small_%d.csv', uint32(posixtime(datetime())));
+writetable(small_table, name)
+name = sprintf('conservatism_large_%d.csv', uint32(posixtime(datetime())));
+writetable(large_table, name)

graphs/circular_graph.m

@@ -0,0 +1,43 @@
+function G = circular_graph(nvert, nedge, directed)
+%CIRCULAR_GRAPH Generates a Matlab graph object that represents a circular
+%graph with a given numer ob vertices.
+%   This function genenerates a Matlab graph object that contains a
+%   circular graph with the given amount of vertices. It can be configured
+%   if the connections should be directed or not and to how many of the
+%   next vertices the connection should be established.
+%
+%   Arguments:
+%       nvert    -> Number of vertices
+%       nedge    -> Number of forward edges per vertex
+%       directed -> [optional] Directed edges or not
+
+if nargin <= 1 || isempty(nedge)
+    nedge = 1;
+elseif nedge >= nvert
+    error('So many neighbours are not possible!')
+end
+if nargin <= 2
+    directed = true;
+end
+
+% Define closure that keeps the index in the ring [1, nvert]
+ring = @(i) mod(i-1, nvert) + 1;
+
+%% Generate graph structure
+A = zeros(nvert);
+for i = 1:nvert
+    A(i, ring(i+(1:nedge))) = 1;
+    
+    % If not direction, also add the other direction
+    if ~directed
+        A(i, ring(i-(1:nedge))) = 1;
+    end
+end
+
+% G needs only to be a digraph if it is directed
+if directed
+    G = digraph(A);
+else
+    G = graph(A);
+end
+end

graphs/laplace_matrix.m

@@ -0,0 +1,9 @@
+function L = laplace_matrix(G)
+%LAPLACE_MATRIX Calculate the Laplacian of the given graph
+%   This function is a replacement for the Matlab built-in laplacian()
+%   function that also works for digraphs, unlike the Matlab version.
+
+A = adjacency(G);
+D = diag(sum(A,2));
+L = D - A;
+end

graphs/triangle_graph.m

@@ -0,0 +1,48 @@
+function G = triangle_graph(baselength)
+%TRIANGLE_GRAPH Generates a Matlab graph object that represents a
+%triangular graph.
+%   This function generates a Matlab graph object that represents a
+%   triangular graph. The generated graphs are always undirected and can be
+%   configured only by their size, which is specified by giving the number
+%   of vertices in the bottommost row.
+
+%% Setup intermediate variables
+% Calculate total number of agents
+N = baselength * (baselength+1) / 2;
+
+% Initial distributing algorithm for the first row
+row_start = 1; % Index of the first node of this row
+row_end   = 1; % Index of the last node of this row
+row_no    = 1; % Number of the row as counted from the top
+
+%% Main algorithm
+A = zeros(N);
+
+for i = 1:N
+    % Check if node is on the bottom row, and if not, add two edges
+    if row_no ~= baselength
+        A(i, i+row_no)   = 1;
+        A(i, i+row_no+1) = 1;
+    end
+    
+    % Check if node is on the left edge of the triangle, if not, add two
+    % edges
+    if i ~= row_start
+        A(i, i-row_no) = 1;
+        A(i, i-1)      = 1;
+    end
+    
+    % Check if node is on the right edge. If so, switch to the next row,
+    % otherwise add two more edges
+    if i ~= row_end
+        A(i, i-row_no+1) = 1;
+        A(i, i+1)        = 1;
+    else
+        row_no    = row_no  + 1;
+        row_start = row_end + 1;
+        row_end   = row_end + row_no;
+    end
+end
+
+G = graph(A);
+end

performance/performance_decomposed.m

@@ -0,0 +1,100 @@
+function [H2, Q, solver_stats] = performance_decomposed(sysD, sysC, L0, p)
+%PERFORMANCE_DECOMPOSED Calculate an upper bound on the H2-norm of a
+%decomposable jump system in a scalable manner
+%   This function calculates an upper bound on the H2-norm of a
+%   decomposable jump system in a scalable manner by applying a decoupled
+%   analysis approach that is described in Theorem 7 of the accompanying
+%   paper.
+%
+%   Arguments:
+%       sysD -> Decoupled part of the system
+%       sysC -> Coupled part of the system
+%       L0   -> Laplacian describing the nominal graph of the system
+%       p    -> Probability of a successful package transmission
+%   Returns:
+%       H2           -> Upper bound on the H2-norm of the system
+%       Q            -> Storage function matrix of the solution
+%       solver_stats -> Timing information about the algorithm
+
+% Setup the SDP solver
+% The offset is used to convert strict LMI into nonstrict ones by pushing
+% the solution away from the boundary.
+offset = 1e-8;
+opts   = sdpsettings('verbose', 0);
+
+% Allocate storage for time measurements
+solver_stats        = struct;
+solver_stats.prep   = NaN;
+solver_stats.solver = NaN;
+solver_stats.total  = NaN;
+solver_stats.yalmip = NaN;
+
+% Initialize return values with reasonable defaults
+H2    = inf;
+timer = tic;
+
+%% Prepare the system description
+[Ac, Bc, Cc, Dc] = ssdata(sysC);
+[Ad, Bd, Cd, Dd] = ssdata(sysD);
+
+nx = size(Ac,1);
+nu = size(Bc,2);
+N  = size(L0,1);
+
+lambda = eig(L0);
+
+%% Solve the SDP from Theorem 7
+Q    = sdpvar(nx, nx);
+Z    = sdpvar(nu, nu, N-1);
+cost = 0;
+Constraints = Q >= offset * eye(nx);
+
+% The eigenvalues of L0 are sorted by Matlab. By iteration only over 2:N,
+% we ignore lambda_1 = 0 and thus the uncontrollable and marginally stable
+% modal subsystem.
+for i = 2:N
+    li  = lambda(i);
+    lit = p*(1-p)*li;
+    Zi  = Z(:,:,i-1);
+    
+    Acl = Ad + p*li*Ac;
+    Bcl = Bd + p*li*Bc;
+    Ccl = Cd + p*li*Cc;
+    Dcl = Dd + p*li*Dc;
+    
+    LMI = Acl'*Q*Acl + Ccl'*Ccl - Q  + 2*lit * (Ac'*Q*Ac + Cc'*Cc);
+    TRC = Bcl'*Q*Bcl + Dcl'*Dcl - Zi + 2*lit * (Bc'*Q*Bc + Dc'*Dc);
+    cost = cost + trace(Zi);
+    
+    % For certain LMIs that are obviously infeasible, Yalmip will refuse to
+    % construct the SDP and issue an error. This try catch will convert
+    % that error into a warning so that we can successfully finish running
+    % the remainder of the script.
+    try
+        Constraints = [ Constraints                     ,...
+                        LMI <= -offset * eye(size(LMI)) ,...
+                        TRC <= -offset * eye(size(TRC)) ];
+    catch ME
+        warning(ME.message)
+        Q = [];
+        return
+    end
+end
+       
+solver_stats.prep = toc(timer);
+sol = optimize(Constraints, cost, opts);
+
+if sol.problem ~= 0
+    warning('YALMIP return an error: %s', sol.info)
+    Q = [];
+else
+    % We calculate gamma^2 with the LMI constraints, so we need to take the
+    % square root here.
+    H2 = sqrt(value(cost));
+    
+    Q  = value(Q);
+    solver_stats.yalmip = sol.yalmiptime;
+    solver_stats.solver = sol.solvertime;
+    solver_stats.total  = toc(timer);
+end
+end