%      SyncGMSC.M
%   
%  Figures to demonstrate the synchronisation problem
%  cf. SyncDemo.m in NCSAufl2
%  for 2nd edition: 12.6.2021
%
% J. Lunze
% 
% 10.9.2023
% Version v. 24.11.2023
% 
%
echo off
clear
close all

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  Harmonic oscillator
figure(1)
%   New oscillator model
A=[0 3; -3 0];
B=[2; 2];      
C=[1 10];
D=0;
n=2;
Oscillator=ss(A, B, C, D);
rlocprintblue(Oscillator);
latexxlabel('Re');
latexylabel('Im');
rotateY;
axis([-4 4 -4 4])
axis('square');
plot([-6 6],[0 0],'k:', [0 0],[-5 5],'k:');
latextitle('$${\tt SyncGMSC:}$$ Root locus of the stable oscillator');
% rlocfind(Oscillator)    %  stable for k=0 ... 0.334
Eig=rlocus(Oscillator,0.167);
plot(real(Eig), imag(Eig), 'kd');
latextext(0.2, 0.4, '$$k_{\rm crit}=0.167$$');
ks=0.1;
Eig=rlocus(Oscillator,ks);
plot(real(Eig), imag(Eig), 'kp');
latextext(-3, 2, '$$k=0.1$$');
xticks([-4 -2 0 2 4]);
yticks([-4 -2 0 2 4]);
hold off

%%   Leader-follower structure and leaderless structure
%
%
figure(2)
N1=4;
L1=[0 0 0 0; -1 1 0 0; 0 -1 1 0; 0 0 -1 1];
Abar1=kron(diag(ones(N1,1)), A)-kron(ks*L1, B*C);
Bbar=zeros(N1*n, 1);
Cbar=kron(diag(ones(N1,1)), C);
Dbar=zeros(N1,1);
OscillatorNetwork=ss(Abar1, Bbar, Cbar, Dbar);
% initial states  
x0=[1; 1; zeros(2*(N1-1),1)];
Tend=10;
dt=0.01;
Time=0:dt:Tend;
[Y1, Time]=initial(OscillatorNetwork, x0, Time);
subplot(4,1,[1 2])
plot(Time, Y1, 'b', Time, Y1(:,1), 'w--');
hold on
plot([0 Tend],[0 0],'k:');
axis([0 Tend -19 19]);
%latexxlabel('$$t$$');
latexylabel('$$y_i$$');
rotateY;
latextitle('$${\tt SyncGMSC:}$$ 4 harmonic oscillators');
xticksempty([0 5 10 15 20]);
xleer;
yticks([-10 0 10]);
hold off
% %   leaderless structure
% L2=[1 0 0 -1; -1 1 0 0; 0 -1 1 0; 0 0 -1 1];
% ks2=ks/2;    % because the largest eigenvalue is 2
% Abar2=kron(diag(ones(N1,1)), A)-kron(ks2*L2, B*C);
% OscillatorNetwork2=ss(Abar2, Bbar, Cbar, Dbar); 
% [Y2, Time]=initial(OscillatorNetwork2, x0, Time);
% subplot(6,1,[3 4])
% plot(Time, Y2, 'b');
% hold on
% plot([0 Tend],[0 0],'k:');
% axis([0 Tend -16 16]);
% xleer;
% latexylabel('$$y_i$$');
% rotateY;
% xticksempty([0 5 10 15 20]);
% yticks([-10 0 10]);
% hold off
%   leaderless structure
L3=[2 -1 0 -1; -1 2 -1 0; 0 -1 2 -1; -1 0 -1 2];   % weight-balanced
ks3=ks/4;    % because the largest eigenvalue is 4
Abar3=kron(diag(ones(N1,1)), A)-kron(ks3*L3, B*C);
OscillatorNetwork3=ss(Abar3, Bbar, Cbar, Dbar); 
[Y3, Time]=initial(OscillatorNetwork3, x0, Time);
subplot(4,1,[3 4])
plot(Time, Y3, 'b', Time, 0.25*sum(Y3'), 'b');
hold on
plot(Time, 0.25*sum(Y3'), 'w--');
plot([0 Tend],[0 0],'k:');
axis([0 Tend -19 19]);
latexxlabel('$$t$$');
%xleer;
latexylabel('$$y_i$$');
rotateY;
xticks([0 5 10 15 20]);
yticks([-10 0 10]);
hold off
%
%  normalised left eigenvector
what=ones(1,N1)*inv(L3+ones(N1,N1));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Convergence of the consensus dynamics
%
figure(3)
N1=9;
AG=[0 1 0 0 0 0 0 1 0 ;
    1 0 0 0 0 1 0 0 0 ;
    0 0 0 1 0 0 1 0 0 ;
    0 0 1 0 1 0 0 1 0 ;
    0 0 0 1 0 0 0 0 0 ;
    0 1 0 0 0 0 0 1 0 ;
    0 0 1 0 0 0 0 0 1 ;
    1 0 0 1 0 1 0 0 0 ;
    0 0 0 0 0 0 1 0 0 ];
L=diag(sum(AG))-AG;
[V,D]=eig(L);
[ds,ind] = sort(diag(D));
Vs = V(:,ind);
lambda2=ds(2);
v2=Vs(:,2);
%  initial state for which the slowest convergence rate is obtained
alpha=1;
x0=alpha*v2;
ACons=-L;
BCons=zeros(N1,1);
CCons=zeros(1,N1);
ConsSyst=ss(ACons, BCons, CCons, 0);
TendCons=5/lambda2;
[Y,TimeCons,XCons]=initial(ConsSyst,x0,TendCons);
subplot(4,1,[1,2]);
plot(TimeCons,XCons,'b');
latexylabel('$$y_i\;\;$$');
rotateY;
xleer;
axis([0 25 -0.6 0.6]);
yticks([-0.5 0 0.5]);
xticksempty([0 5 10 15 20 25 30]);
latextitle('$${\tt SyncGMSC:}$$ Convergence of consensus dynamics');
%  another initial state
%SyncGMSC2=rng;                         % save s to restore the generator setting
%save('SyncGMSC2', 'SyncGMSC2');
load('SyncGMSC2');       % Seed (initialisation) of the number generator
rng(SyncGMSC2);
x02=x0+rand(N1,1);
[Y,TimeCons2,XCons2]=initial(ConsSyst,x02,TendCons);
subplot(4,1,[3,4]);
plot(TimeCons2,XCons2,'b');
latexylabel('$$y_i\;\;$$');
rotateY;
latexxlabel('$$t$$');
axis([0 25 -0.6 1.2]);
yticks([-0.5 0 0.5 1 1.5]);
xticks([0 5 10 15 20 25 30]);

%%  Convergence of consensus systems with random graph
%   cf. RandomGraphDemo.m 
%   The convergence is evaluated by all states reaching a bound of
%   epsilonbar around the average consensus value
%
figure(4)
N1=20;
Tend=20;
dt=0.01;
Time=0:dt:Tend;
x0=[N1; zeros(N1-1,1)];
xbar=1;
kmax=200;
p1=0.20;
p2=0.35;
p3=0.8;
epsilon=0.001;    % bound for lambda2 to indicate connectivity of the graph
epsilonbar=0.05;  % tolerance of convergence
pbar=log(N1)/(N1-1);   % lower bound for connected random graphs
%Prob=[0.02:0.01:1];
Prob=[0.17:0.005:0.4 0.45:0.05:1];
kp=0;
ProbUsed=[];
Lambda2=[];       % average algebraic connectivity
                  % in dependence upon the probability p
Convergence=[];   % Convergence time (T_5%)
Timeconst=[];
for p=Prob
    if p>pbar
        p
        lambda2Sum=0;
        ConvSum=0;
        k2=0;
        %  Determine for the given p the algebraic connectivity lambda2
        %  provided that the graph is connected
        %  use kmax realisations
        for k1=1:kmax
            %  Generation of the adjacency matrix of a random graph
            A=ceil(rand(N1,N1)-(1-p)*ones(N1,N1)); 
            %  make the matrix symmetric for undirected graphs
            for k3=1:N1
                for k4=k3+1:N1
                    A(k3,k4)=A(k4,k3);
                end
            end
            A=A-diag(diag(A));     % set the main diagonal elements zero
            %  Generate the Laplacian matrix
            D=diag(sum(A));      % degree matrix
            L=D-A;               % Laplacian matrix
            %  determine the number of vanishing eigenvalues
            Eig=eig(L);
            EigSorted=sort(Eig,1,'ascend');
            if EigSorted(2)>epsilon      % if the graph is connected 
                %  determine the time instant at which the bound is reached
                ConsSyst=ss(-L, zeros(N1,1), zeros(1,N1), 0);
                [Y,Time,X]=initial(ConsSyst,x0,Time);
                for k3=length(Time):-1:1
                    if max(abs(X(k3,:)-xbar))<epsilonbar
                       ConvTime=Time(k3);
                    end
                end
                ConvSum=ConvSum+ConvTime;
                k2=k2+1;
                lambda2Sum=lambda2Sum+EigSorted(2);
            end
        end
        kp=kp+1;
        Lambda2(kp)=lambda2Sum/k2;
        ProbUsed(kp)=p;
        Convergence(kp)=ConvSum/k2;
        NumberOfSamples(kp)=k2;  % number of connected graphs obtained
        if abs(p-p1)<epsilon
            Time1=ConvSum/k2;
        end
        if abs(p-p2)<epsilon
            Time2=ConvSum/k2;
        end
        if abs(p-p3)<epsilon
            Time3=ConvSum/k2;
        end
    end
end
%
subplot(4,1,[1,2]);
plot(ProbUsed, Lambda2, 'b');
hold on
plot([pbar, pbar],[0 N1],'b:');
plot(ProbUsed, ProbUsed*N1, 'b--');
latextext(0.5, 13.5, '$$pN$$');
latexylabel('$$\bar{\lambda}_2\;\;$$');
rotateY;
xleer;
axis([0 1 0 N1]);
yticks([0 10 20 30]);
xticksempty([0 0.2 0.4 0.6 0.8 1]);
%latextext(0.8, 3.4, ['$$N=\;$$', num2str(N1)]); 
latextitle('$${\tt SyncGMSC:}$$ Settling time of the consensus dynamics');
hold off
%
subplot(4,1,[3,4]);
plot(ProbUsed, Convergence, 'b');
hold on
plot([pbar, pbar],[0 N1],'b:');
plot(p1,Time1,'bo',...
            'MarkerEdgeColor','b',...
            'MarkerFaceColor','b');
plot(p2,Time2,'bo',...
            'MarkerEdgeColor','b',...
            'MarkerFaceColor','b');
plot(p3,Time3,'bo',...
            'MarkerEdgeColor','b',...
            'MarkerFaceColor','b');
latexylabel('$$\bar{T}_{5 \%}\;$$');
rotateY;
latexxlabel('$$p$$');
latextext(0.8*pbar,-1.1,'$$\bar{p}$$');
axis([0 1 0 9]);
yticks([0 3 6 9]);
xticks([0 0.2 0.4 0.6 0.8 1]);
% for k1=1:length(Lambda2)
%     Timeconst(k1)=1/Lambda2(k1);  
% end
% plot(ProbUsed, 4*Timeconst, 'b--');  % precisely the same as Convergence
hold off

%%  Demonstration of the result by three examples
%
figure(5)
N1=20;
Tend=6;     % simulation time horizon
dt=0.01;
Time=0:dt:Tend;
%SyncGMSC1=rng;                         % save s to restore the generator setting
%save('SyncGMSC1', 'SyncGMSC1');
load('SyncGMSC1');       % Seed (initialisation) of the number generator
rng(SyncGMSC1);
for p=[p1,p2,p3]
%  Generation of the adjacency matrix of a random graph
            A=ceil(rand(N1,N1)-(1-p)*ones(N1,N1)); 
            %  make the matrix symmetric for undirected graphs
            for k3=1:N1
                for k4=k3+1:N1
                    A(k3,k4)=A(k4,k3);
                end
            end
            A=A-diag(diag(A));     % set the main diagonal elements zero
            %  Generate the Laplacian matrix
            D=diag(sum(A));      % degree matrix
            L=D-A;               % Laplacian matrix
            %  determine the number of vanishing eigenvalues
            Eig=eig(L);
            EigSorted=sort(Eig,1,'ascend');
            if EigSorted(2)>epsilon      % if the graph is connected 
                %  determine the consensus dynamics
                Tend=6;     % simulation time horizon
                ConsSyst=ss(-L, zeros(N1,1), zeros(1,N1), 0);
                [Y,Time,X]=initial(ConsSyst,x0,Time);
                if p==p1
                    subplot(6,1,[1,2]);
                    k3=length(X);
                    while max(abs(X(k3,:)-xbar))<epsilonbar
                        k3=k3-1;
                    end
                    Teps1=k3*dt;
                    A1=A;
                end
                if p==p2
                    subplot(6,1,[3 4]);
                    k3=length(X);
                    while max(abs(X(k3,:)-xbar))<epsilonbar
                        k3=k3-1;
                    end
                    Teps2=k3*dt;
                    A2=A;
                end
                if p==p3
                    subplot(6,1,[5 6]);
                    k3=length(X);
                    while max(abs(X(k3,:)-xbar))<epsilonbar
                        k3=k3-1;
                    end
                    Teps3=k3*dt;
                    A3=A;
                end
                plot(Time,X,'b');
                hold on
            end
end
subplot(6,1,[1,2]);
plot([Teps1 Tend],[1.05 1.05],'r', [Teps1 Tend],[0.95 0.95],'r');
plot([Teps1 Teps1],[2 0],'r');
latextext(Teps1,5,'$$T_{5 \%}$$');
latexylabel('$$x_i\;$$');
rotateY;
xleer;
yticks([0 10 20 30]);
xticksempty([0 1 2 3 4 5 6]);
axis([0 Tend 0 22]);
latextext(3,17,'$$p_1=0.20$$');
hold off
%
subplot(6,1,[3, 4]);
plot([Teps2 Tend],[1.05 1.05],'r', [Teps2 Tend],[0.95 0.95],'r');
plot([Teps2 Teps2],[2 0],'r');
latextext(Teps2,5,'$$T_{5 \%}$$');
latexylabel('$$x_i\;$$');
rotateY;
xleer;
yticks([0 10 20 30]);
xticksempty([0 1 2 3 4 5 6]);
axis([0 Tend 0 22]);
latextext(3,17,'$$p_2=0.35$$');
hold off
%
subplot(6,1,[5 6]);
plot([Teps3 Tend],[1.05 1.05],'r', [Teps3 Tend],[0.95 0.95],'r');
plot([Teps3 Teps3],[2 0],'r');
latextext(Teps3,5,'$$T_{5 \%}$$');
latexylabel('$$x_i\;$$');
rotateY;
latexxlabel('$$t$$');
yticks([0 10 20 30]);
xticks([0 1 2 3 4 5 6]);
axis([0 Tend 0 22]);
latextext(3,17,'$$p_3=0.80$$');
hold off



%%  Convergence of consensus systems with random graph
%   use the normalised Laplacian matrix
%
figure(6)
N1=20;
Tend=20;
dt=0.01;
Time=0:dt:Tend;
x0=[N1; zeros(N1-1,1)];
xbar=1;
kmax=200;
p1=0.20;
p2=0.35;
p3=0.8;
epsilon=0.001;    % bound for lambda2 to indicate connectivity of the graph
epsilonbar=0.05;  % tolerance of convergence
pbar=log(N1)/(N1-1);   % lower bound for connected random graphs
%Prob=[0.02:0.01:1];
Prob=[0.17:0.005:0.4 0.45:0.05:1];
kp=0;
ProbUsed=[];
Lambda2=[];       % average algebraic connectivity
                  % in dependence upon the probability p
Convergence=[];   % Convergence time (T_5%)
Timeconst=[];
for p=Prob
    if p>pbar
        lambda2Sum=0;
        ConvSum=0;
        k2=0;
        %  Determine for the given p the algebraic connectivity lambda2
        %  provided that the graph is connected
        %  use kmax realisations
        for k1=1:kmax
            %  Generation of the adjacency matrix of a random graph
            A=ceil(rand(N1,N1)-(1-p)*ones(N1,N1)); 
            %  make the matrix symmetric for undirected graphs
            for k3=1:N1
                for k4=k3+1:N1
                    A(k3,k4)=A(k4,k3);
                end
            end
            A=A-diag(diag(A));     % set the main diagonal elements zero
            %  Determine the normalised adjacency matrix
            for k3=1:N1
                if sum(A(k3,:))~=0
                    Ahat(k3,:)=A(k3,:)/sum(A(k3,:));
                else
                    Ahat(k3,:)=A(k3,:);
                end
            end
            %  Generate the Laplacian matrix
            Dhat=diag(sum(Ahat'));      % degree matrix
            Lhat=Dhat-Ahat;               % Laplacian matrix
            %  determine the number of vanishing eigenvalues
            Eig=eig(Lhat);
            EigSorted=sort(Eig,1,'ascend');
            if EigSorted(2)>epsilon      % if the graph is connected 
                %  determine the time instant at which the bound is reached
                ConsSyst=ss(-Lhat, zeros(N1,1), zeros(1,N1), 0);
                [Y,Time,X]=initial(ConsSyst,x0,Time);
                xbar=X(end,1);
                for k3=length(Time):-1:1
                    if max(abs(X(k3,:)-xbar))<epsilonbar
                       ConvTime=Time(k3);
                    end
                end
                ConvSum=ConvSum+ConvTime;
                k2=k2+1;
                lambda2Sum=lambda2Sum+EigSorted(2);
            end
        end
        kp=kp+1;
        Lambda2(kp)=lambda2Sum/k2;
        ProbUsed(kp)=p;
        Convergence(kp)=ConvSum/k2;
        NumberOfSamples(kp)=k2;  % number of connected graphs obtained
        if abs(p-p1)<epsilon
            Time1=ConvSum/k2;
        end
        if abs(p-p2)<epsilon
            Time2=ConvSum/k2;
        end
        if abs(p-p3)<epsilon
            Time3=ConvSum/k2;
        end
    end
end
%
subplot(4,1,[1,2]);
plot(ProbUsed, Lambda2, 'b');
hold on
plot([pbar, pbar],[0 N1],'b:');
plot(ProbUsed, ProbUsed, 'b--');
latexylabel('$$\bar{\lambda}_2\;\;$$');
rotateY;
xleer;
axis([0 1 0 1]);
yticks([0 0.5 1]);
xticksempty([0 0.2 0.4 0.6 0.8 1]);
%latextext(0.8, 3.4, ['$$N=\;$$', num2str(N1)]); 
latextitle('$${\tt SyncGMSC:}$$ Settling time of the consensus dynamics');
latextext(0.2, 0.8, '$$\hat{L}$$');
hold off
%
subplot(4,1,[3,4]);
plot(ProbUsed, Convergence, 'b');
hold on
plot([pbar, pbar],[0 N1],'b:');
plot(p1,Time1,'bo',...
            'MarkerEdgeColor','b',...
            'MarkerFaceColor','b');
plot(p2,Time2,'bo',...
            'MarkerEdgeColor','b',...
            'MarkerFaceColor','b');
plot(p3,Time3,'bo',...
            'MarkerEdgeColor','b',...
            'MarkerFaceColor','b');
latexylabel('$$\bar{T}_{5 \%}\;$$');
rotateY;
latexxlabel('$$p$$');
latextext(0.8*pbar,-2.1,'$$\bar{p}$$');
axis([0 1 0 15]);
yticks([0 5 10 15 20 25]);
xticks([0 0.2 0.4 0.6 0.8 1]);
% for k1=1:length(Lambda2)
%     Timeconst(k1)=1/Lambda2(k1);  
% end
% plot(ProbUsed, 4*Timeconst, 'b--');  % precisely the same as Convergence
hold off



%%  Demonstration of the result by three examples
%
figure(7)
N1=20;
Tend=25;     % simulation time horizon
dt=0.01;
Time=0:dt:Tend;
%SyncGMSC1=rng;                  % save s to restore the generator setting
%save('SyncGMSC1', 'SyncGMSC1');
load('SyncGMSC1');       % Seed (initialisation) of the number generator
rng(SyncGMSC1);
for p=[p1,p2,p3]
%  Generation of the adjacency matrix of a random graph
            A=ceil(rand(N1,N1)-(1-p)*ones(N1,N1)); 
            %  make the matrix symmetric for undirected graphs
            for k3=1:N1
                for k4=k3+1:N1
                    A(k3,k4)=A(k4,k3);
                end
            end
            A=A-diag(diag(A));     % set the main diagonal elements zero
            %  Determine the normalised adjacency matrix
            for k3=1:N1
                if sum(A(k3,:))~=0
                    Ahat(k3,:)=A(k3,:)/sum(A(k3,:));
                else
                    Ahat(k3,:)=A(k3,:);
                end
            end
            %  Generate the Laplacian matrix
            Dhat=diag(sum(Ahat'));      % degree matrix
            Lhat=Dhat-Ahat;               % Laplacian matrix
            %  determine the number of vanishing eigenvalues
            Eig=eig(Lhat);
            EigSorted=sort(Eig,1,'ascend');
            if EigSorted(2)>epsilon      % if the graph is connected 
                %  determine the consensus dynamics
                ConsSyst=ss(-Lhat, zeros(N1,1), zeros(1,N1), 0);
                [Y,Time,X]=initial(ConsSyst,x0,Time);
                if p==p1
                    subplot(4,1,[1,2]);
                    xbar1=X(end,1);           %  consensus value
                    k3=length(X);
                    while max(abs(X(k3,:)-xbar1))<epsilonbar
                        k3=k3-1;
                    end
                    Teps1=k3*dt;
                else
                    subplot(4,1,[3 4]);
                    xbar2=X(end,1);           %  consensus value
                    k3=length(X);
                    while max(abs(X(k3,:)-xbar2))<epsilonbar
                        k3=k3-1;
                    end
                    Teps2=k3*dt;
                end
                plot(Time,X,'b');
                hold on
            end
end
subplot(4,1,[1,2]);
plot([Teps1 Tend],[1.05 1.05]*xbar1,'r', [Teps1 Tend],[0.95 0.95]*xbar1,'r');
plot([Teps1 Teps1],[2 0],'r');
latextext(Teps1,-1.8,'$$T_{5 \%}$$');
latexylabel('$$x_i\;$$');
rotateY;
xleer;
yticks([0 10 20 30]);
xticksempty([0 5 10 15 20 25 30]);
axis([0 Tend 0 20]);
latextext(3,17,'$$p=0.20$$');
hold off
%
subplot(4,1,[3 4]);
plot([Teps2 Tend],[1.05, 1.05]*xbar2 ,'r', [Teps2 Tend],[0.95 0.95]*xbar2,'r');
plot([Teps2 Teps2],[2 0],'r');
latextext(Teps2,3.3,'$$\bar{T}_{5 \%}$$');
latexylabel('$$x_i\;$$');
rotateY;
latexxlabel('$$t$$');
yticks([0 10 20 30]);
xticks([0 5 10 15 20 25 30]);
axis([0 Tend 0 20]);
latextext(3,17,'$$p=0.35$$');
hold off

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  Multirotor fleet
%
figure(8)
radius=100;   % 100 m
N1=40;
%x0=10*(1-2*rand(N1,1));
x0 =[
    5.9032
    7.2280
   -7.3334
   -5.7000
   -6.6451
   -4.0865
   -1.3817
   -9.9485
   -1.0263
   -3.2857
    7.3393
    0.4447
    7.6851
    2.9674
   -4.4065
   -6.9372
   -3.4741
   -7.0646
   -5.4798
    3.3420
   -9.7386
   -5.8189
    8.1235
    2.2640
    4.4023
    6.6874
   -1.1659
   -5.5218
   -4.5198
   -9.4069
    7.0121
   -9.8173
   -2.9809
    4.5212
   -6.4076
    3.1193
    0.2474
   -1.6539
    3.3830
    6.3698];
Tend=250;
dt=0.1;
Time=0:dt:Tend;
kmax=100;
epsilonbar=0.05;  % tolerance of convergence
Prob=[0:0.005:0.2 0.22:0.02:0.3];
p1=0.05;
kp=0;
ProbUsed=[];
Convergence=[];   % Convergence time (T_5%)
for p=Prob
        p
        ConvSum=0;
        k2=0;
        %  use kmax realisations
        for k1=1:kmax
            %  Generation of the adjacency matrix of a random graph
            A=ceil(rand(N1,N1)-(1-p)*ones(N1,N1)); 
            %  make the matrix symmetric for undirected graphs
            for k3=1:N1
                for k4=k3+1:N1
                    A(k3,k4)=A(k4,k3);
                end
            end
            A=A-diag(diag(A));     % set the main diagonal elements zero
            % add the neighbouring couplings
            for k3=1:N1-1
                A(k3+1,k3)=1;
                A(k3,k3+1)=1;
            end
            A(1,N1)=1;
            A(N1,1)=1;
             %  Generate the Laplacian matrix
            D=diag(sum(A'));      % degree matrix
            L=D-A;                % Laplacian matrix
            %  determine the time instant at which the bound is reached
            ConsSyst=ss(-L, zeros(N1,1), zeros(1,N1), 0);
            [Y,Time,X]=initial(ConsSyst,x0,Time);
            xbar=X(end,1);
            for k3=length(Time):-1:1
                if max(abs(X(k3,:)-xbar))<epsilonbar
                   ConvTime=Time(k3);
                end
            end
            k2=k2+1;
            ConvSum=ConvSum+ConvTime;
        end
        kp=kp+1;
        ProbUsed(kp)=p;
        Convergence(kp)=ConvSum/k2;
        if abs(p-p1)<0.001
            Time1=ConvSum/k2;
        end
end
%
subplot(4,1,[1,2]);
plot(ProbUsed, Convergence, 'b');
hold on
plot(p1,Time1,'bo',...
            'MarkerEdgeColor','b',...
            'MarkerFaceColor','b');
latexylabel('$$\bar{T}_{5 \%}\;$$');
rotateY;
latexxlabel('$$p$$');
axis([0 0.3 0 80]);
yticks([0 20 40 60 80]);
xticks([0 0.05 0.1 0.15 0.2 0.25 0.3]);
hold off


%%  Demonstration of the result by two examples
%
figure(9)
p0=0;
N1=40;
Tend=150;
dt=0.5;     % for dt>0.2: eps-figure impossible
Time=0:dt:Tend;
%SyncGMSC3=rng;                         % save s to restore the generator setting
%save('SyncGMSC3', 'SyncGMSC3');
load('SyncGMSC3');       % Seed (initialisation) of the number generator
rng(SyncGMSC3);
for p=[p0,p1]
 %  Generation of the adjacency matrix of a random graph
            A=ceil(rand(N1,N1)-(1-p)*ones(N1,N1)); 
            %  make the matrix symmetric for undirected graphs
            for k3=1:N1
                for k4=k3+1:N1
                    A(k3,k4)=A(k4,k3);
                end
            end
            A=A-diag(diag(A));     % set the main diagonal elements zero
            % add the neighbouring couplings
            for k3=1:N1-1
                A(k3+1,k3)=1;
                A(k3,k3+1)=1;
            end
            A(1,N1)=1;
            A(N1,1)=1;
             %  Generate the Laplacian matrix
            D=diag(sum(A'));      % degree matrix
            L=D-A;                % Laplacian matrix
            %  determine the time instant at which the bound is reached
            ConsSyst=ss(-L, zeros(N1,1), zeros(1,N1), 0);
            [Y,Time,X]=initial(ConsSyst,x0,Time);
            if p==p0
                    subplot(4,1,[1 2]);
                    xbar0=X(end,1);
                    k3=length(X);
                    while max(abs(X(k3,:)-xbar0))<epsilonbar
                        k3=k3-1;
                    end
                    Teps0=k3*dt;
                    A0=A;
                end
            if p==p1
                    subplot(4,1,[3 4]);
                    xbar1=X(end,1);
                    k3=length(X);
                    while max(abs(X(k3,:)-xbar1))<epsilonbar
                        k3=k3-1;
                    end
                    Teps1=k3*dt;
                    A1=A;
             end
             plot(Time,X+radius,'b');
             hold on
             plot([0 Tend],[radius radius],'k--');
end
subplot(4,1,[1 2]);
plot([Teps0 Tend],radius+[1.05 1.05]*xbar0,'r', [Teps0 Tend],radius+[0.95 0.95]*xbar0,'r');
plot([Teps0 Teps0],radius+[2 0]*xbar0,'r');
latextext(Teps0,radius+4,'$$T_{5 \%}$$');
latexylabel('$$x_i$$ in m');
rotateY;
xleer;
yticks([90 100 110 120]);
xticksempty([0 50 100 150 200]);
axis([0 150 90 110]);
latextext(100,108,'$$p_0=0$$');
hold off
%
subplot(4,1,[3 4]);
plot([Teps1 Tend],radius+[1.05 1.05]*xbar1,'r', [Teps1 Tend],radius+[0.95 0.95]*xbar1,'r');
plot([Teps1 Teps1],radius+[2 0]*xbar1,'r');
latextext(Teps1,radius+4,'$$T_{5 \%}$$');
latexylabel('$$x_i$$ in m');
rotateY;
latexxlabel('$$t$$ in s');
yticks([90 100 110 120]);
xticks([0 50 100 150 200]);
axis([0 150 90 110]);
latextext(100,108,'$$p_1=0.05$$');
hold off

%%  Graphical representation of the two communication graphs
%
figure(12);
G0=graph(A0);
G1=graph(A1);
subplot(2,3,[1 2])
%axis('square');
%subplot(3,2,1)
G.Nodes.Name=[];
h=plot(G0,'Layout','circle',...
    'NodeLabel',G.Nodes.Name,...
    'MarkerSize',5,'LineWidth',1.6);
xleer;
yleer;
%
subplot(2,3,[4 5])
%subplot(3,2,2)
axis('square');
G.Nodes.Name=[];
h=plot(G1,'Layout','circle',...
    'NodeLabel',G.Nodes.Name,...
    'MarkerSize',5,'LineWidth',1.6);
xleer;
yleer;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%   Synchronisation of three oscillators
%  Bidirectional couplings by the resistor R
%   cf. LCNetworkDemo.m  in NCSAufl2
%
figure(10)
N1=3;
%   Oscillators
L1=0.1;   % 100 mH
C1=0.00001;   % 10 muF
A=[0 1/C1;  -1/L1 0];
B1=[-1/C1; 0];
C1=[-1  0];
D=0;
n=2;
Oscillator=ss(A, B1, C1, D);
VC01=-1;  % 1 V
IL01=0;  % 0 A
x01=[VC01; IL01];
Tend=0.015;
dt=0.0001;
Time=[0:dt:Tend];
%
kc1=0.005;    % hence, R1=1/kc = 10 Ohm
kc2=0.04;
L=[1 -1  0;
   -1  2 -1;
   0 -1  1 ]; 
%
subplot(4,1,[1,2]);
L2=L*kc1; 
Abar=kron(diag(ones(N1,1)), A)-kron(L2, B1*C1);
Bbar=zeros(N1*n, 1);
Cbar=kron(diag(ones(N1,1)), C1);
Network=ss(Abar, Bbar, Cbar, 0);
x0=[x01; zeros((N1-1)*n, 1)];
Y=initial(Network, x0, Time);
plot(Time, Y, 'b');
latexylabel('$$V_i$$ in V');
rotateY;
latexxlabel('$$t$$ in s');
latextitle('$${\tt SyncGMSC:}$$ Oscillator network behaviour');
axis([0 Tend -1.1 1.1]);
yticks([-1 0 1]);
xleer
xticksempty([0 0.005 0.01 0.015]);
latextext(0.01, 0.8, ['$$1/R=$$' num2str(kc1)]);
%
subplot(4,1,[3,4]);
L3=L*kc2;  % hence, R1=1/kc2 = 1.25 Ohm
Abar=kron(diag(ones(N1,1)), A)-kron(L3, B1*C1);
Network=ss(Abar, Bbar, Cbar, 0);
Y=initial(Network, x0, Time);
plot(Time, Y, 'b');
latexylabel('$$V_i$$ in V');
rotateY;
latexxlabel('$$t$$ in s');
axis([0 Tend -1.1 1.1]);
yticks([-1 0 1]);
xticks([0 0.005 0.01 0.015]);
latextext(0.01, 0.8, ['$$1/R=$$' num2str(kc2)]);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%   Comparison of the convergence of consensus dynamics
%  a) complete graph:  lambda2=N
%  b) ring graph:      lambda2=0.4679
%  By selecting appropriate feedback gains, the convergence rate becomes
%  the same
%
figure(11)
N1=9;
L1=N1*eye(N1,N1)-ones(N1,N1);
L2=[2 -1 0 0 0 0 0 0 -1;
    -1 2 -1 0 0 0 0 0 0;
    0 -1 2 -1 0 0 0 0 0;
    0 0 -1 2 -1 0 0 0 0;
    0 0 0 -1 2 -1 0 0 0;
    0 0 0 0 -1 2 -1 0 0;
    0 0 0 0 0 -1 2 -1 0;
    0 0 0 0 0 0 -1 2 -1;
    -1 0 0 0 0 0 0 -1 2;];
kc1=1/N1;
kc2=1/0.4679;
x0=[9; 0; 0; 0; 0; 0; 0; 0; 0];
ACons=-kc1*L1;
BCons=zeros(N1,1);
CCons=zeros(1,N1);
ConsSyst=ss(ACons, BCons, CCons, 0);
TendCons=5;
[Y,TimeCons,XCons]=initial(ConsSyst,x0,TendCons);
subplot(4,1,[1,2]);
plot(TimeCons,XCons,'b', 5, 1,'bo');
latexylabel('$$y_i\;\;$$');
rotateY;
xleer;
axis([0 5 -1 10]);
yticks([0 5 10]);
xticksempty([0 1 2 3 4 5]);
latextitle('$${\tt SyncGMSC:}$$ Convergence of consensus dynamics');

%
ACons=-kc2*L2;
ConsSyst2=ss(ACons, BCons, CCons, 0);
[Y,TimeCons2,XCons2]=initial(ConsSyst2,x0,TendCons);
subplot(4,1,[3,4]);
plot(TimeCons2,XCons2,'b', 5, 1,'bo');
latexylabel('$$y_i\;\;$$');
rotateY;
latexxlabel('$$t$$');
axis([0 5 -1 10]);
yticks([0 5 10]);
xticks([0 1 2 3 4 5]);


%%  figure(12) used


%%   Figures
%
epsfigc15('SyncGMSC');
