%      RandUQNet.M
%   Störverhalten von Teilsysteme mit zufällig
%   gewählten Kopplungen
%   from   RandUQNetzStat.m  of 20.3.2010
%
% J. Lunze
% 
% 22.2.2010
% Version vom 28.6.2010
% 
%
echo off
clear
close all
kmax=300;               % Anzahl der untersuchten Netzwerke für dieselbe Wahrscheinlichkeit
P=[0.005:0.01:0.6, 1];    % Verbindungswahrscheinlichkeit
XT=1;                   % Transformatorreaktanz
N1=200;                  % Anzahl der Knoten

%% Analyse eines Widerstandsnetzwerks
%  Statische Analyse des ungeregelten Netzes
%
figure(1)
subplot(3,1,[1,2]);
%  ungeregelte Kraftwerke
%  untersucht wird die Situation, in der am Knoten 1 die Stromeinspeisung
%  so groß gewählt wird, dass die Spannung U_1=1 entsteht
%  Dafür wird untersucht, an wie vielen anderen Knoten die
%  Spannungsänderung eine Schranke (0.05, 0.1, 0.2, 0.3) überschreitet
%
%  Admittanzmatrix Adjbest enthält die Transformatorreaktanz 1/XT
%
%  zufällige Kopplungen zwischen den Knoten sind mit der Wahrscheinlichkeit
%  p vorhanden
%
%
%%%%%%%%%%%%%%%%%%
kN=0;
for p=P
    p                %  Kontrolle des Berechnungsfortschritts
    Nd1=0; Nd2=0; Nd3=0; Nd4=0;
for k=1:kmax
% 
%  Die Koppelmatrix entsteht aus zufällig gewählten Verbindungen. 
%  Es wird angenommen, dass jede Verbindung den Leitwert 1 hat.
ZufM=ceil(rand(N1+1,N1+1)-(1-p)*ones(N1+1,N1+1)); % Erzeugen einer zuf. Matrix
                                % die folgenden Schritte machen die Matrix
                                % symmetrisch:
for k1=1:N1+1
    for k2=k1+1:N1+1
        ZufM(k1, k2)=ZufM(k2, k1);
    end
end
ZufM=ZufM-diag(diag(ZufM));     % Nullsetzen der Hauptdiagonalelemente
Summe=sum(ZufM);                % spaltenweises Summieren
Adj=-ZufM+diag(Summe);          % Hauptdiagonalelemente ist Summe der
                                % anderen Spaltenelemente 
                                % alle anderen Elemente mit negativem
                                % Vorzeichen
Adjbest=Adj(1:N1, 1:N1);          % Streichen des Bezugsknotens: 
                                % bestimmte Admittanzmatrix Adjbest
                                % Hinzufügen der Ableitung aller Knoten 
                                % zum Bezugsknoten mit Leitwert 1/XT
Adjbest=Adjbest+1/XT*diag(ones(1,N1));
%
%  Berechnung der Knotenspannungen U_2 ... U_N
%
U=-inv(Adjbest(2:N1, 2:N1))*Adjbest(2:N1, 1);
%  Nd = Anzahl der Teilsysteme, deren Spannung u_i
%       größer als ubar=0.05, 0.1, 0.2, 0.3 ist
Nd1=Nd1+sum(ceil(U-0.05*ones(N1-1,1)));
Nd2=Nd2+sum(ceil(U-0.1*ones(N1-1,1)));
Nd3=Nd3+sum(ceil(U-0.2*ones(N1-1,1)));
Nd4=Nd4+sum(ceil(U-0.3*ones(N1-1,1)));
end
kN=kN+1;
Ndmittel1(kN)=Nd1/kmax;
Ndmittel2(kN)=Nd2/kmax;
Ndmittel3(kN)=Nd3/kmax;
Ndmittel4(kN)=Nd4/kmax;
end
semilogx(P, Ndmittel1, 'b');
hold on
semilogx(P, Ndmittel2, 'b');
semilogx(P, Ndmittel3, 'b');
%semilogx(P, Ndmittel4, 'b-.');
plot([0.13 0.13],[200 0] ,'b--');
latextext(0.13, -15, '$$\bar{p}=0.13$$');
axis([0 1 0 N1+1]);
latexxlabel('$$p$$');
latexylabel('$$n$$');
rotateY;
xticks([0.001 0.01 0.1 1]);
yticks([0 50 100 150 200])
latextitle('$${\tt RandUQNet:}$$ Mean number of disturbed vertices');
hold off

%%  Ersatzwiderstand eines Netzwerkes mit zufälliger Struktur
%
figure(2)
subplot(3,1,[1,2]);
kN=0;
for p=P
    p                %  Kontrolle des Berechnungsfortschritts
    RersSum=0;
    k2=0;
    for k=1:kmax
        %  from RandomGraphDemo.m:
        %  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
        l22=L(2,2);
        l12=L(1,2);
        l1T=L(1,3:N1);
        l2T=L(2,3:N1);
        l1=L(3:N1,1);
        l2=L(3:N1,2);
        L22=L(3:N1,3:N1);
        if det(L22)~=0
           k2=k2+1;
           RersSum=RersSum+1/(l22-l2T*inv(L22)*l2);
        end
    end
kN=kN+1;
Rers(kN)=RersSum/k2;
end
plot(P, Rers,'b');
latexxlabel('$$p$$');
latexylabel('$$R_{\rm eq}$$');
rotateY;
xticks([0 0.2 0.4 0.6 0.8 1]);
yticks([0 0.1 0.2 0.3 0.4])
axis([0 1 0 0.4]);
latextitle('$${\tt RandUQNet:}$$ Resistance of a random network');
hold off

%%     Bilder auf .EPS- und .PS-Dateien
%
epsfigc('RandUQNet');



