%   Script NetwAnal.m
%
%   Electrical network analysis
%
% J. Lunze
% 6.10.2022
% Version of 14.10.2022
%
%   
clear
close all


%%  Fundamental cycles - Exercise
%
C=[1 0 0
   -1 0 0
   0 0 1
   1 0 1
   0 0 -1
   -1 -1 0
   0 1 0
   0 1 0];
N=[-1 -1 0 0 0 0 0 0;
    0 0 0 -1 -1 -1 -1 0
    0 1 -1 1 0 0 0 0
    1 0 0 0 0 1 0 1
    0 0 0 0 0 0 1 -1
    0 0 1 0 1 0 0 0];



%%  Branch current analysis
%
%figure(1)
N1=4;
M1=6;
C=[1 0 0; 0 1 -1; -1 1 0; -1 0 1; 0 -1 0; 0 0 1];
N=[1 0 1 0 1 0; 0 1 -1 1 0 0; 0 -1 0 0 -1 -1; -1 0 0 -1 0 1];
Ntilde=N(1:N1-1,:);
inv([C'; Ntilde])



%% Node voltage analysis
%
Lperp=[3 -1 -1; -1 3 -1; -1 -1 3];


%%  Kron reduction
%
N1=6;
L=-[0 1 0 1 0 0 0;
    1 0 1 0 1 0 0
    0 1 0 0 0 1 0
    1 0 0 0 1 0 1
    0 0 1 1 0 1 1
    0 0 1 0 1 0 1
    0 0 0 1 1 1 0];
L=L-diag(sum(L'));
Lperp=L(1:N1-1, 1:N1-1);
L12=L(1:N1-1,N1);
L21=L(N1,1:N1-1);
L22=L(N1,N1);
Lred=Lperp-L12*inv(L22)*L21


%%   Example 9.7
%
G1=1;
G2=1;
G3=1;
G4=1;
G5=1;
L=[G1+G2 0 -G2 -G1; 0 G4+G5 -G5 -G4; -G2 -G5 G2+G3+G5 -G3; -G1 -G4 -G3 G1+G3+G4];
% delete the vertices 3 and 4 
L11=L(1:2,1:2);
L12=L(1:2,3:4);
L21=L(3:4,1:2);
L22=L(3:4,3:4);
Lred=L11-L12*inv(L22)*L21;

%%  Kron reduction of a latter
%
AG=[0 0 1 0 0 0 0 0 0 0;
    0 0 0 1 0 0 0 0 0 0;
    1 0 0 1 1 0 0 0 0 0;
    0 1 1 0 0 1 0 0 0 0;
    0 0 1 0 0 1 1 0 0 0;
    0 0 0 1 1 0 0 1 0 0;
    0 0 0 0 1 0 0 1 1 0;
    0 0 0 0 0 1 1 0 0 1;
    0 0 0 0 0 0 1 0 0 1;
    0 0 0 0 0 0 0 1 1 0];
L=diag(sum(AG'))-AG;
% delete the vertices 9 and 10
L11=L(1:8,1:8);
L12=L(1:8,9:10);
L21=L(9:10,1:8);
L22=L(9:10,9:10);
Lred=L11-L12*inv(L22)*L21;
%
% delete the vertices 7, 8, 9, 10
L11=L(1:6,1:6);
L12=L(1:6,7:10);
L21=L(7:10,1:6);
L22=L(7:10,7:10);
Lred=L11-L12*inv(L22)*L21;


%%  Exercise: Kron reduction of a ladder
%
AG=[0 0 1 0 0 0 0 0 0 0 0 0;
    0 0 0 0 0 1 0 0 0 0 0 0;
    1 0 0 1 0 0 1 0 0 0 0 0;
    0 0 1 0 1 0 0 1 0 0 0 0;
    0 0 0 1 0 1 0 0 1 0 0 0;
    0 1 0 0 1 0 0 0 0 1 0 0;
    0 0 1 0 0 0 0 1 0 0 1 0;
    0 0 0 1 0 0 1 0 1 0 1 0;
    0 0 0 0 1 0 0 1 0 1 0 1;
    0 0 0 0 0 1 0 0 1 0 0 1;
    0 0 0 0 0 0 1 1 0 0 0 1;
    0 0 0 0 0 0 0 0 1 1 1 0];
L=diag(sum(AG'))-AG;
% delete the vertices 11 and 12
L11=L(1:10,1:10);
L12=L(1:10,11:12);
L21=L(11:12,1:10);
L22=L(11:12,11:12);
Lred=L11-L12*inv(L22)*L21;
%   effective resistance
P=10;
Lplus=inv(Lred+1/P*ones(P,P))-1/P*ones(P,P);
R12=Lplus(1,1)+Lplus(2,2)-2*Lplus(1,2);
%
%  delete the vertices 3, 4, 5, 6
% permutation matrix to get the 
% vertex order: 1, 2, 7, 8, 9, 10, 11, 12; 3, 4, 5, 6
P=[1 0 0 0 0 0 0 0 0 0 0 0;
   0 1 0 0 0 0 0 0 0 0 0 0;
   0 0 0 0 0 0 0 0 1 0 0 0;
   0 0 0 0 0 0 0 0 0 1 0 0;
   0 0 0 0 0 0 0 0 0 0 1 0;
   0 0 0 0 0 0 0 0 0 0 0 1;
   0 0 1 0 0 0 0 0 0 0 0 0;
   0 0 0 1 0 0 0 0 0 0 0 0;
   0 0 0 0 1 0 0 0 0 0 0 0;
   0 0 0 0 0 1 0 0 0 0 0 0;
   0 0 0 0 0 0 1 0 0 0 0 0;
   0 0 0 0 0 0 0 1 0 0 0 0];
Lneu=P'*L*P;
% delete the last four vertices
L11=Lneu(1:8,1:8);
L12=Lneu(1:8,9:12);
L21=Lneu(9:12,1:8);
L22=Lneu(9:12,9:12);
Lred=L11-L12*inv(L22)*L21;
%  effective resistance
P=8;
Lplus=inv(Lred+1/P*ones(P,P))-1/P*ones(P,P);
R12=Lplus(1,1)+Lplus(2,2)-2*Lplus(1,2);

%%   Exercise: Kron reduction
%
figure(1)
R2=1;
R3=1/2;
R4=1/2;
R5=1/3;
R6=1/3;
%
%   Triangle-star transformation
RA=(R2*R4+R2*R5+R4*R5)/R2;
RB=(R2*R4+R2*R5+R4*R5)/R5;
RC=(R2*R4+R2*R5+R4*R5)/R4;
RB3=RB*R3/(RB+R3);
RC6=RC*R6/(RC+R6);
Reff=RA*(RB3+RC6)/(RA+RB3+RC6);
%
Reffvar=[];
R5var=0.1:0.1:6;
k1=0;
for R5n=R5var
    k1=k1+1;
    RA=(R2*R4+R2*R5n+R4*R5n)/R2;
    RB=(R2*R4+R2*R5n+R4*R5n)/R5n;
    RC=(R2*R4+R2*R5n+R4*R5n)/R4;
    RB3=RB*R3/(RB+R3);
    RC6=RC*R6/(RC+R6);
    Reffvar(k1)=RA*(RB3+RC6)/(RA+RB3+RC6);
end
subplot(2,1,1)
plot(R5var, Reffvar, 'b');
latexxlabel('$$R_5$$');
latexylabel('$$R_{\rm eff}$$');
rotateY;
yticks([0.2 0.4 0.6]);
xticks([0 2 4 6 8 10]);
%
%  Kron reduction
G2=1/R2;
G3=1/R3;
G4=1/R4;
G5=1/R5;
G6=1/R6;
L=[G3+G4   0     -G4     -G3;
     0   G5+G6   -G5     -G6
    -G4   -G5  G2+G4+G5  -G2
    -G3   -G6    -G2   G2+G3+G6];
Lred=L(1:2,1:2)-L(1:2,3:4)*inv(L(3:4,3:4))*L(3:4,1:2);
Rred=1/Lred(1,1);



%%   Exercise: Kron reduction
%
AG=[0 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0;
    1 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0;
    1 0 0 1 0 0 1 0 0 0 0 0 0 0 0 0;
    1 0 1 0 1 0 0 1 0 0 0 0 0 0 0 0;
    0 1 0 1 0 1 0 0 1 0 0 0 0 0 0 0;
    0 1 0 0 1 0 0 0 0 1 0 0 0 0 0 0;
    0 0 1 0 0 0 0 1 0 0 1 0 0 0 0 0;
    0 0 0 1 0 0 1 0 0 0 0 1 0 0 0 0;
    0 0 0 0 1 0 0 0 0 1 0 0 1 0 0 0;
    0 0 0 0 0 1 0 0 1 0 0 0 0 1 0 0;
    0 0 0 0 0 0 1 0 0 0 0 1 0 0 1 0;
    0 0 0 0 0 0 0 1 0 0 1 0 1 0 1 0;
    0 0 0 0 0 0 0 0 1 0 0 1 0 1 0 1;
    0 0 0 0 0 0 0 0 0 1 0 0 1 0 0 1;
    0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 1;
    0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 0];
L=diag(sum(AG'))-AG;
% permutation matrix to get the 
% vertex order: 8, 9, 1, 2,..., 7, 10, ..., 16
P=[0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0;
   0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0;
   0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0;
   0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0;
   0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0;
   0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0;
   0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0;
   1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
   0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
   0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0;
   0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0;
   0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0;
   0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0;
   0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0;
   0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0;
   0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1;];
Lneu=P'*L*P;
% determine the effective resistance with respect to the first two vertices
L11=Lneu(1:2,1:2);
L12=Lneu(1:2,3:16);
L21=Lneu(3:16,1:2);
L22=Lneu(3:16,3:16);
Lred=L11-L12*inv(L22)*L21;
Reff=1/Lred(1,1);
%
%  delete the edge (14 - 16)
AG=[0 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0;
    1 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0;
    1 0 0 1 0 0 1 0 0 0 0 0 0 0 0 0;
    1 0 1 0 1 0 0 1 0 0 0 0 0 0 0 0;
    0 1 0 1 0 1 0 0 1 0 0 0 0 0 0 0;
    0 1 0 0 1 0 0 0 0 1 0 0 0 0 0 0;
    0 0 1 0 0 0 0 1 0 0 1 0 0 0 0 0;
    0 0 0 1 0 0 1 0 0 0 0 1 0 0 0 0;
    0 0 0 0 1 0 0 0 0 1 0 0 1 0 0 0;
    0 0 0 0 0 1 0 0 1 0 0 0 0 1 0 0;
    0 0 0 0 0 0 1 0 0 0 0 1 0 0 1 0;
    0 0 0 0 0 0 0 1 0 0 1 0 1 0 1 0;
    0 0 0 0 0 0 0 0 1 0 0 1 0 1 0 1;
    0 0 0 0 0 0 0 0 0 1 0 0 1 0 0 0;
    0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 1;
    0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0]
%  repeat the determination of Reff
L=diag(sum(AG'))-AG;
P=[0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0;
   0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0;
   0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0;
   0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0;
   0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0;
   0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0;
   0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0;
   1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
   0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
   0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0;
   0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0;
   0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0;
   0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0;
   0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0;
   0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0;
   0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1;];
Lneu=P'*L*P;
L11=Lneu(1:2,1:2);
L12=Lneu(1:2,3:16);
L21=Lneu(3:16,1:2);
L22=Lneu(3:16,3:16);
Lred=L11-L12*inv(L22)*L21;
Reff2=1/Lred(1,1);
%  delete the vertex 6 and three edges
AG=[0 1 1 1 0  0 0 0 0 0 0 0 0 0 0;
    1 0 0 0 1  0 0 0 0 0 0 0 0 0 0;
    1 0 0 1 0  1 0 0 0 0 0 0 0 0 0;
    1 0 1 0 1  0 1 0 0 0 0 0 0 0 0;
    0 1 0 1 0  0 0 1 0 0 0 0 0 0 0;
%    0 1 0 0 1  0 0 0 1 0 0 0 0 0 0;
    0 0 1 0 0  0 1 0 0 1 0 0 0 0 0;
    0 0 0 1 0  1 0 0 0 0 1 0 0 0 0;
    0 0 0 0 1  0 0 0 1 0 0 1 0 0 0;
    0 0 0 0 0  0 0 1 0 0 0 0 1 0 0;
    0 0 0 0 0  1 0 0 0 0 1 0 0 1 0;
    0 0 0 0 0  0 1 0 0 1 0 1 0 1 0;
    0 0 0 0 0  0 0 1 0 0 1 0 1 0 1;
    0 0 0 0 0  0 0 0 1 0 0 1 0 0 0;
    0 0 0 0 0  0 0 0 0 1 1 0 0 0 1;
    0 0 0 0 0  0 0 0 0 0 0 1 0 1 0];
%  repeat the determination of Reff
L=diag(sum(AG'))-AG;
P=[0 0 1 0 0  0 0 0 0 0 0 0 0 0 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 0 0 0 0 0 0;
   0 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 0 0 0 0;
 %  0 0 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 0 0 0;
   1 0 0 0 0  0 0 0 0 0 0 0 0 0 0;
   0 1 0 0 0  0 0 0 0 0 0 0 0 0 0;
   0 0 0 0 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 0;
   0 0 0 0 0  0 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 0 0 0 0  0 0 0 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 0 0  0 0 0 0 0 0 0 0 0 1;];
Lneu=P'*L*P;
L11=Lneu(1:2,1:2);
L12=Lneu(1:2,3:15);
L21=Lneu(3:15,1:2);
L22=Lneu(3:15,3:15);
Lred=L11-L12*inv(L22)*L21;
Reff3=1/Lred(1,1);


%%
epsfigc('NetwAnal');
    