%      RandomNetDemo.M
%   
%  Consensus in undirected random graphs
%  [Hatano/Mesbahi IEEE Trans. AC 2005]
%
%   probabilities p >0.1
%   only a few switches until convergence
%   continuous-time simulation
%
% J. Lunze
% 
% 11.1.2018
% Version v. 27.11.2018
% for 2nd edition: 16.4.2021
%
echo off
clear
close all
P=[0.04:0.02:1, 1];      % Probability of edges
%P=[0.001:0.002:0.1];    % ... for details with small p
kmax=300;                % Number of graphs investigated for the same p

%%   Consensus in random graphs
%
%   xdot = -L x           first-order integrator agents
%   xbar = average of x0
%
%   Laplacian matrix changes every time step
%   Ts - sampling time
%   x(k+1) = P(k) x(k)    with P(k) = expm(-L(k)*Ts)
figure(1)
N1=12;                  % number of agents
p1=0.08;
p=p1;
Ts=1;
ke=8;
dt=0.01;                 % step size for continuous-time simulation
kstep=Ts/dt;
Time=[0:dt:ke*Ts];
X=[];                    % continuous-time state
Xd=[];                   % state at sampling times
V=[];                    % Lyapunov function
Vd=[];                   % Lyapunov function at sampling times
No=[];                   % number of communication links
%RandomNet=rng;                         % save s to restore the generator setting
%save('RandomNet', 'RandomNet');
load('RandomNet');       % Seed (initialisation) of the number generator
rng(RandomNet);
%
%   Multi-agent system
%x0=abs(rand(N1,1));   % random initial state
%x0=x0/sum(x0)*N1;        %  nonnegative initial state in [0, 1]
x0=[1.1593
    1.4511
    1.3621
    1.7787
    0.3655
    0.7770
    0.8191
    0.2092
    1.1034
    1.0904
    0.1776
    1.7066];
xbar=sum(x0)/N1;      %  consensus state
X(1,:)=x0;
Xd(1,:)=x0;
c=sum(x0);
V(1)=x0'*x0-c^2/N1;   % Lyapunov function
Vd(1)=V(1);
ksim=0;                 % simulation time index
Lambda2=[];             % lamda_2(P)
for k1=1:ke
    %  Generation of the adjacency matrix of a random graph
    G=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
            G(k3, k4)=G(k4, k3);
        end
    end
    G=G-diag(diag(G));     % set the main diagonal elements zero
    %  Generate the Laplacian matrix
    D=diag(sum(G));      % degree matrix
    No(k1)=sum(sum(D))/2;  % Number of communication links used
    L=D-G;               % Laplacian matrix
    %
    %   Continuous-time simulation
    Ac=expm(-L*dt);
    for k2=1:kstep
        ksim=ksim+1;
        xnew=Ac*X(ksim,:)';     % next state
        X(ksim+1,:)=xnew;
        V(ksim+1)=xnew'*xnew-c^2/N1;
    end
    %
    %   Discrete-time simulation
    Ad=expm(-L*Ts);
    xdnew=Ad*Xd(k1,:)';     % next state
    Xd(k1+1,:)=xdnew;
    Vd(k1+1)=xdnew'*xdnew-c^2/N1;
    Lambda=eig(Ad);
    Lambda=sort(Lambda,1,'descend');
    Lambda2(k1)=Lambda(2);
end
subplot(3,1,[1,2]);
plot(Time, X,'b');
hold on
plot([0:ke]*Ts, Xd, 'bo');
plot([0 ke]*Ts,[xbar, xbar],'k:');
latextitle(['$${\tt RandomNetDemo:}$$ Consensus behaviour for $$p=$$' num2str(p)]);
%latexxlabel('$$k$$');
xleer;
latexylabel('$$x_i$$');
rotateY;
yticks([0 1 2]);
axis([0 ke 0 2]);
hold off
%
%  Lyapunov function
subplot(3,1,3)
plot(Time, V,'b');
hold on
plot([0:ke]*Ts, Vd, 'bo');
latexxlabel('$$t$$');
latexylabel('$$V$$');
rotateY;
xticks([0 1 2 3 4 5 6 7 8 9 ]);
yticks([0 2 4]);
axis([0 ke 0 4]);
hold off
%
%   Connectivity
figure(2)
mbar=N1*(N1+1)/2*p;
subplot(2,1,1)
stairs([0:ke], [No No(end)],'LineWidth',1.5,'Color', 'b');
%plot([0:ke-1]+0.5, No, 'b');  % the number of edges holds true for the
                              % time period between consecutive samplings
hold on
plot([0 ke],[mbar mbar],'k--');
latexxlabel('$$t$$');
latexylabel('$$m\;$$');
rotateY;
axis([0 ke 0 15]);
yticks([0 5 10 15]);
xticks([0 1 2 3 4 5 6 7 8 9 10]);
latextitle(['$${\tt RandomNetDemo:}$$ Consensus behaviour for $$p=$$' num2str(p)]);
hold off
% subplot(2,1,2)
% plot([1:ke]*Ts, Lambda2,'k');
% hold on;
% plot([1:ke]*Ts, Lambda2,'ko');
% latexxlabel('$$t$$');
% latexylabel('$$V$$');
% rotateY;
% xticks([20 30 40]);
% yticks([1 2]);
% axis([0 ke 0.5 2]);
% hold off

%%   Analysis of time-to-consensus
%
figure(3)
N1=20;
P=[0.04:0.02:0.5, 0.52:0.05:1, 1];
epsilon=0.001;    % Precision of consensus
%
for qk=[0.3 0.5 0.8]
kProb=0;
Ttc=[];
for p=P
    p        % show iteration of p
    kProb=kProb+1;
    kend=0;
for kAverage=1:kmax     % Average over kmax graphs
    X=[];    
%   Multi-agent system
x0=abs(rand(N1,1));   
x0=x0/sum(x0)*N1;        %  nonnegative initial state in [0, 1]
X(1,:)=x0;
V2=max(x0)-min(x0);
k1=0;
while V2>epsilon
    k1=k1+1;
    %  Generation of the adjacency matrix of a random graph
    G=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
            G(k3, k4)=G(k4, k3);
        end
    end
    G=G-diag(diag(G));     % set the main diagonal elements zero
    %  Generate the normed Laplacian matrix
    Ahat=inv(diag(max(sum(G'),ones(1,N1))))*G; % normalised adjacency matrix
    Lhat=diag(sum(Ahat'))-Ahat;                % normalised Laplacian matrix
    Ad=eye(N1,N1)-qk*Lhat;
    xnew=Ad*X(k1,:)';
    X(k1+1,:)=xnew;
    V2=max(xnew)-min(xnew);
end
kend=kend+k1;
end
Ttc(kProb)=kend/kmax;      % time to consensus
end
subplot(3,1,[1,2])
plot(P, Ttc, 'b');
hold on;
end
latextitle(['RandomNetDemo: Time-to-consensus $$\epsilon=$$' num2str(epsilon)]);
latexxlabel('$$p$$');
latexylabel('$$k_\varepsilon \;$$');
rotateY;
xticks([0 0.2 0.4 0.6 0.8 1]);
yticks([0 10 20 30 40 50]);
axis([0 1 0 max(Ttc*Ts)]);
latextext(0.6, 30, '$$q=0.3$$');
latextext(0.68, 15, '$$0.5$$');
%latextext(0.68, 10, '$$0.8$$');
hold off



%%  Random agreement
%   use the graphs of "switching system consensus" in ConsensusDemo.m
%   figure(20)
%
figure(6)
N1=5;
A1=[0 0 0 0 1; 0 0 0 0 0; 0 0 0 0 0; 0 0 0 0 0; 1 0 0 0 0];
A2=[0 0 0 1 0; 0 0 0 0 0; 0 0 0 0 1; 1 0 0 0 0; 0 0 1 0 0];
A3=[0 1 1 0 0; 1 0 1 0 0; 1 1 0 0 0; 0 0 0 0 0; 0 0 0 0 0];
D1=diag(sum(A1));
D2=diag(sum(A2));
D3=diag(sum(A3));
L1=D1-A1;
L2=D2-A2;
L3=D3-A3;
Ts=0.05;
P1=expm(-L1*Ts);
P2=expm(-L2*Ts);
P3=expm(-L3*Ts);
%
p1=0.2;       % probability of using the graph G1
p2=0.5;
p3=1-p1-p2;
kend=90;
Time=[0:kend-1];
Prob=rand(kend,1);
%
x0=[0; 1; 0; 1; 0];
xbar=sum(x0)/N1;
X=[];
X(1,:)=x0';
Index=[];
V=[];         % Lypunov function
V(1,1)=x0'*x0-N1*xbar*xbar;
for k1=2:kend
    if Prob(k1)<=p1
        P=P1;
        Index(k1-1)=1;
    elseif Prob(k1)<=p1+p2
        P=P2;
        Index(k1-1)=2;
    else
        P=P3;
        Index(k1-1)=3;
    end
     X(k1,:)=P*X(k1-1,:)';
     V(k1,1)=X(k1,:)*X(k1,:)'-N1*xbar*xbar;
end
%
subplot(4,1,1)
plot([0:kend-2], Index, 'bo');
latexylabel('$$i\;$$');
rotateY;
xleer;
axis([0 kend+1 0.9 3.1]);
yticks([1 2 3]);
xticksempty([0 20 40 60 80 100]);
%
subplot(4,1,[2,3])
stairs(Time,X,'Color','b','LineWidth',1.5);
hold on
xleer;
latexylabel('$$x_i$$');
rotateY;
axis([0 kend+1 -0.05 1.05]);
yticks([0 0.5 1]);
xticksempty([0 20 40 60 80 100]);
%  consensus value
xbar=sum(x0)/N1;
plot(kend+1, xbar, 'bo');
latextext(kend+4, xbar, '$$\bar{x}$$');
hold off;
%
%  Convergence analysis
ExpP=p1*P1+p2*P2+p3*P3;
eig(ExpP);
%
% Lyapunov function
subplot(4,1,4)
stairs(Time,V,'Color','b','LineWidth',1.5);
hold on
%plot(Time,V,'bo');
latexxlabel('$$k$$');
latexylabel('$$V$$');
rotateY;
axis([0 kend+1 0 1.2]);
yticks([0 1]);
xticks([0 20 40 60 80 100]);
hold off;



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Random synchronisation of oscillators
%
figure(7)
%   oscillator model from SyncDemo.m
A=[0 3; -3 0];
B=[1; 1];      
C=[1 10];
D=0;
n=2;
%  parameters
N1=12;                  % number of agents
Cbar=kron(eye(N1,N1), C);
p=0.1;                  % connection probability of the random graph
kP=0.01;
Ts=0.1;                   % sampling time
dt=0.1;                 % simulation time
kstep=ceil(Ts/dt);
Tend=20;
kend=ceil(Tend/Ts);
TimeS=[0:dt:kend*Ts];
X=[];                    % continuous-time state
Y=[];
RandomNet2=rng;                         % save s to restore the generator setting
save('RandomNet2', 'RandomNet2');
load('RandomNet2');       % Seed (initialisation) of the number generator
rng(RandomNet2);
x0=abs(rand(N1*n,1));   % random initial state
x0=x0/sum(x0)*N1;        %  nonnegative initial state in [0, 1]
X(1,:)=x0;
Y(1,:)=Cbar*x0;
ksim=0;
for k1=1:kend
    %  Generation of the adjacency matrix of a random graph
    G=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
            G(k3, k4)=G(k4, k3);
        end
    end
    G=G-diag(diag(G));     % set the main diagonal elements zero
    %  Generate the Laplacian matrix
    D=diag(sum(G));      % degree matrix
    L=D-G;               % Laplacian matrix
    Abar=kron(eye(N1,N1), A)- kP*kron(L, B*C);
    %
    %   Continuous-time simulation
    Ac=expm(-Abar*dt);
    for k2=1:kstep
        ksim=ksim+1;
        xnew=Ac*X(ksim,:)';     % next state
        X(ksim+1,:)=xnew;
     end
end
subplot(3,1,[1,2]);
plot(TimeS, X,'b');
hold on
latextitle(['$${\tt RandomNetDemo:}$$ Synchronisation behaviour for $$p=$$'...
                            num2str(p), '$$, \;T_{\rm s}=\;$$' num2str(Ts)]);
%latexxlabel('$$k$$');
xleer;
latexylabel('$$y_i$$');
rotateY;
yticks([-1 0 1]);
axis([0 Tend -1.3 1.3]);
hold off



%%   Figures
epsfigc15('RandomNetDemo');

