diff --git a/CheckAndersLMI.m b/CheckAndersLMI.m
new file mode 100644
index 0000000000000000000000000000000000000000..c14096169553d1e4288a956f0c6a897eeaa1af11
--- /dev/null
+++ b/CheckAndersLMI.m
@@ -0,0 +1,195 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% This code simulates distributed minimax adaptive control algorithm for
+% uncertain networked systems.
+%
+% Copyrights Authors: 1) Venkatraman Renganathan - Lund University, Sweden.
+% 2) Anders Rantzer - Lund University, Sweden.
+% 3) Olle Kjellqvist - Lund University, Sweden.
+%
+% Email: venkatraman.renganathan@control.lth.se
+% anders.rantzer@control.lth.se
+% olle.kjellqvist@control.lth.se
+%
+% Date last updated: 20 October, 2023.
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+% Make a fresh start
+clear all; close all; clc;
+
+% set properties for plotting
+set(groot,'defaultAxesTickLabelInterpreter','latex');
+set(groot,'defaulttextinterpreter','latex');
+set(groot,'defaultLegendInterpreter','latex');
+
+% Flag deciding whether to generate new data or load existing data
+% dataPrepFlag = 1: Generates new network data
+% dataPrepFlag = 0: Loads existing network data
+dataPrepFlag = 1;
+
+% When dataPrepFlag = 1: Generate new data for network
+if(dataPrepFlag)
+
+ disp('Generating new data for network as dataPrepFlag = 1.');
+
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ % Define the network data using graph structure
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+ % Specify Number of nodes
+ N = 6;
+ % Specify Number of edges
+ M = N - 1;
+ % Create a random graph with N nodes and M edges
+ [adjMatrix, I, graphVariable] = GenerateRandomTreeGraph(N);
+ % Get number of edges on graph
+ numEdges = size(graphVariable.Edges.EndNodes, 1);
+ fprintf('Generated random graph is acyclic and has only one component \n');
+
+ % Plot the graph - large graph can be time consuming & may not be informative
+ if N <= 300
+ plot(graphVariable);
+ else
+ disp('The graph is too large to display.');
+ end
+
+ % Infer state and input dimensions
+ [n,m] = size(I);
+
+ % Get the degree vector
+ degreeVec = sum(adjMatrix, 1)';
+
+ % Get one of the node indeces with maximum degree for plotting
+ [~, maxDegreeNodeIndex] = max(degreeVec);
+
+ % Form input matrix B with parameter b and matrix I such that B/b = I
+ b = 0.1;
+ B = b*I;
+
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ % Prepare the dynamics matrix A for each Model
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+ % Specify the number of possible plant network parameters
+ numModels = 2;
+
+ % Populate all the possible plant models
+ for dIter = 1:numModels
+ % Placeholders for a matrix for each component
+ ai = zeros(n, 1);
+ % Flag for stopping while loop
+ dynamicsPreparationFlag = 1;
+ % Counter for populating ai vector
+ iter = 1;
+ % Loop through until ai^2 + bb^{T} < ai is satisfied
+ while dynamicsPreparationFlag
+ % Generate a random ai
+ ai(iter, 1) = rand;
+ % Check condition
+ if(ai(iter, 1)^2 - ai(iter, 1) + 2*b^2*degreeVec(iter, 1) < 0)
+ iter = iter + 1;
+ end
+ % if counter exceeds limit, form the A matrix & stop while loop
+ if(iter > n)
+ A = diag(ai);
+ dynamicsPreparationFlag = 0;
+ end
+ end
+ % Prepare A and B matrices
+ AMatrices{dIter} = A;
+ BMatrices{dIter} = B;
+
+ end
+
+ AVector = zeros(n, numModels);
+ for dIter = 1:numModels
+ AVector(:, dIter) = diag(AMatrices{dIter});
+ end
+ Arowvectors = {};
+ for i = 1:N
+ Arowvectors{i} = AVector(i, :);
+ end
+ Acombs = combvec(Arowvectors{:}).';
+ % Find min & max a values for each node from all possible numModels values
+ minAVector = min(AVector, [], 2);
+ maxAVector = max(AVector, [], 2);
+ % Form the Max A matrix using maxAVector
+ maxAMatrix = diag(maxAVector);
+ % Find the network level global minimum and maximum a value
+ netMinA = min(minAVector);
+ netMaxA = max(maxAVector);
+
+ % Infer the number of possible plants
+ numPlants = numModels^(N);
+
+ % Compute Hinfty control for each possible plants
+ for i = 1:numPlants
+ APossibles{i} = diag(Acombs(i, :));
+ BPossibles{i} = B;
+ % Compute H_infinity control as u = Kx, with K = B'*(A-I)^{-1}
+ KinfPossibles{i} = BPossibles{i}'*inv(APossibles{i} - eye(n));
+ end
+ disp('Finished Computing H_infinity Control Gains');
+
+ % Check Anders LMI for each possible plant
+ gamma = 100;
+ P_non = (2/(1 - netMaxA))*eye(n);
+ P_non_inv = ((1 - netMaxA)/2)*eye(n);
+ P_wt = inv(P_non_inv - gamma^(-2)*eye(n));
+ Q = eye(n);
+ R = eye(M);
+ happy = 1;
+
+ for k = 1:numPlants
+ Ak = APossibles{k};
+ if(happy == 0)
+ break;
+ end
+ for l = 1:numPlants
+ Al = APossibles{l};
+ if(happy == 0)
+ break;
+ end
+ for p = 1:numPlants
+ if(happy == 0)
+ break;
+ end
+ if(k ~= l && l == p)
+ continue;
+ end
+ Kp = KinfPossibles{p};
+ A_cl_kp = Ak + B*Kp;
+ A_cl_lp = Al + B*Kp;
+ AndersRHS = Q + Kp'*R*Kp + 0.25*(A_cl_kp+A_cl_lp)'*P_wt*(A_cl_kp+A_cl_lp) - 0.25*gamma^(2)*(A_cl_kp-A_cl_lp)'*(A_cl_kp-A_cl_lp);
+ AndersLMI = P_non - AndersRHS;
+ if(all(eig(AndersLMI) >= 0))
+ happy = 1;
+ else
+ happy = 0;
+ end
+ end
+ end
+ end
+
+
+ % Finished Checking Anders LMI
+ if(happy == 1)
+ disp('Choice of Pkl matrices satisfy Anders LMI.');
+ else
+ disp('Choice of Pkl matrices DO NOT satisfy Anders LMI.');
+ end
+
+ % Save the generated network data into mat file.
+ disp('Saving the generated small network data.');
+ %save('smallNetworkDynamicsData.mat');
+else
+ disp('Loading existing data for small network as dataPrepFlag = 0.');
+ %load('smallNetworkDynamicsData.mat');
+end
+
+
+
+