Commit f047e026 authored by Jakob Gabriel's avatar Jakob Gabriel
Browse files

misc.synchronisationWithSpectrumConstraints: implements the synchronisation...

misc.synchronisationWithSpectrumConstraints: implements the synchronisation problem of homogenious networks with spectrum constrains via the LMI-toolbox
parent 37a57d89
function [gain, X, Y, dynamics] = synchronisationWithSpectrumConstraints(A, B, laplacian, NameValue)
% misc.synchronisationWithSpectrumConstraints solves the synchronisation problem for i = 1 , ..., N
% homogenious agents
% xi_t = A xi + B ui
% that communicate via a connected network that is descirbed by the laplacian.
% A feedback gain K = Y / X, ui = K xi, is calculated using LMIs, such that the spectrum of
% dynamics = kron(I, A) - kron(laplacian, B*K)
% only contains one copy of the spectrum of A and the other eigenvalues are inside a region that
% can is described with further input parameters "margin", "discRadius", "sectorAngle".
%
% misc.synchronisationWithSpectrumConstraints(A, B, L, "margin" m) ensures that the eigenvalues are
% left of -m
%
% misc.synchronisationWithSpectrumConstraints(A, B, L, "discRadius" r) ensures that the eigenvalues
% are within a disc of radius r.
%
% misc.synchronisationWithSpectrumConstraints(A, B, L, "sectorAngle" a) ensures that the eigenvalues
% are within a sector of angle a.
%
% The constraints margin, discRadius, sector angle can also be combined.
%
% This function is based on the publication
% Seyboth, Georg S., Wei Ren, and Frank Allgöwer. "Cooperative control of linear multi-agent
% systems via distributed output regulation and transient synchronization." Automatica 68
% (2016): 132-139.
% preprint: https://arxiv.org/pdf/1406.0085.pdf
% and uses the matlab LMI tools.
arguments
A double;
B double;
laplacian double;
NameValue.margin (1, 1) double = 0;
NameValue.discRadius (1, 1) double = inf;
NameValue.sectorAngle (1, 1) double = pi;
end % arguments
n = size(A, 1); % length of state vector xi
p = size(B, 2); % length of input ui
M = []; % M and L are parameters for the LMIs that ensure the constraints.
L = [];
if NameValue.margin ~= 0
assert(NameValue.margin >= 0, "margin must be positive to obtain stability")
M = blkdiag(M, 1);
L = blkdiag(L, 2 * NameValue.margin);
end
if NameValue.discRadius < inf
assert(NameValue.discRadius >= 0, "discRadius must be positive to obtain stability")
M = blkdiag(M, [0, 1; 0, 0]);
L = blkdiag(L, -NameValue.discRadius * eye([2, 2]));
end
if NameValue.sectorAngle < pi
assert(NameValue.sectorAngle >= 0, "sectorAngle must be positive to obtain stability")
M = blkdiag(M, [sin(NameValue.sectorAngle), cos(NameValue.sectorAngle); ...
-cos(NameValue.sectorAngle), sin(NameValue.sectorAngle)]);
L = blkdiag(L, zeros([2, 2]));
end
lambda = eig(laplacian);
% it is sufficient to consider
% 1. the smallest non-zero real eigenvalue,
% 2. the largest real eigenvalue,
% 3. all complex-conjugated eigenvalue-pairs.
assert(all(real(lambda) >= 0), "laplacian must have non-negative eigenvalues");
assert(sum(abs(lambda) <= 10*eps) == 1, ...
"graph must be connected, hence the laplacian must have exactly 1 eigenvalue in the origin 0")
lambdaReal = lambda(imag(lambda) == 0);
lambdaReal = lambdaReal(lambdaReal > 0);
l = unique([min(lambdaReal); max(lambdaReal); lambda(imag(lambda)~=0)]);
assert(isreal(l), "LMI solver not yet implemented for laplacians with imaginary eigenvalues" ...
+ " consider https://de.mathworks.com/help/robust/ug/advanced-topics.html")
% use lmi toolbox
setlmis([])
x = lmivar(1,[n, 1]);
y = lmivar(2,[p, n]);
for lt = 1 : numel(l) % set up one lmi for each value in l
for it = 1 : size(M, 1) % set up block matrix for lmi lt
for jt = 1 : size(M, 2)
lmiterm([lt, it, jt, x], M(it, jt) * A, 1, "s");
lmiterm([lt, it, jt, y], M(it, jt) * (-l(lt)) * B, 1, "s");
lmiterm([lt, it, jt, x], L(it, jt), 1);
end % for jt = 1 : size(M, 2)
end % for it = 1 : size(M, 1)
end % for lt = 1 : numel(l)
lmiterm([lt+1, 1, 1, x], -1, 1); % ensures that X > 0 be demanding -X < 0
% calculate and obtain the resultX and Y
lmis = getlmis;
options = zeros(5, 1);
%options(5) = 1;
[tmin,xfeas] = feasp(lmis, options);
X = dec2mat(lmis,xfeas,x);
Y = dec2mat(lmis,xfeas,y);
% if abs(tmin) < 1e-6
% warning("no feasible result for the LMI was found")
% end
% verify result
for lt = 1 : numel(l)
[lhs, rhs] = showlmi(evallmi(lmis, xfeas), lt);
% check if the lhs of the lmi is the one that we intended to solve
lhsIntended = kron(L, X) ...
+ kron(M, A * X - l(lt) * B * Y) ...
+ kron(M.', (A * X - l(lt)' * B * Y).');
assert(max(abs(lhs - lhsIntended), [], "all") <= 1e-6, "probably the wrong lmi was solved");
assert(all(rhs==0, "all"), "probably the wrong lmi was solved");
end % for it = 1 : numel(l)
gain = Y / X;
dynamics = kron(eye(size(laplacian)), A) - kron(laplacian, B * gain);
end % synchronisationWithSpectrumConstraints()
\ No newline at end of file
function [tests] = testSynchronisationWithSpectrumConstraints()
% result = runtests('unittests.misc.testSynchronisationWithSpectrumConstraints')
tests = functiontests(localfunctions);
end
function testExample(tc)
mySubSystems = {[0, 1; -1, 0], [0, 1; -1, 0], [0, 1; 0, 0]};
A = blkdiag(mySubSystems{:});
B = blkdiag([0; 1], [0; 1], [0; 1]);
adjacency = [0, 0, 0; 2, 0, 1; 0, 3, 0];
laplacian = diag(sum(adjacency, 2)) - adjacency;
margin = 1;
discRadius = 10;
sectorAngle = pi/4;
[gain, X, Y, dynamics] = misc.synchronisationWithSpectrumConstraints(A, B, laplacian, ...
"margin", margin, "discRadius", discRadius, "sectorAngle", sectorAngle);
spectrum = eig(dynamics);
spectrumStable = setdiff(...
round(spectrum, 14), ...
round(eig(A), 14));
tc.verifyTrue(all(real(spectrumStable) < - margin));
tc.verifyTrue(all(abs(spectrumStable) < discRadius));
tc.verifyTrue(all(angle(-spectrumStable) < sectorAngle));
end % testExample
\ No newline at end of file
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment