nsw_cluster_smallworld

 Cluster-small world model to generation rand topo power grids
 start from cluster lattice, instead of physical nodal locations.
 Inputs:
   N - network size;
   degModel - degree distribution model;
   mdeg - average nodal degree;
   p_lrwd - probability of rewiring links;
   alph, bta - transition probability in the Markov chain model
 Outputs:
   A - link admittance matrix;
   Ma - connectivity matrix;
   indx - link indexs;
   Md - distance matrix;
   lambda2 - algebraic connecity, i.e., lambda2(L)
   youtN - actual network size of resulting topology;
   L - the Laplacian matrix

 wzf, Dec 2008

0001 function [A,Ma,indx,Md,lambda2,youtN,L] = nsw_cluster_smallworld(N,degModel,mdeg,p_lrwd,alph,bta)
0002 % Cluster-small world model to generation rand topo power grids
0003 % start from cluster lattice, instead of physical nodal locations.
0004 % Inputs:
0005 %   N - network size;
0006 %   degModel - degree distribution model;
0007 %   mdeg - average nodal degree;
0008 %   p_lrwd - probability of rewiring links;
0009 %   alph, bta - transition probability in the Markov chain model
0010 % Outputs:
0011 %   A - link admittance matrix;
0012 %   Ma - connectivity matrix;
0013 %   indx - link indexs;
0014 %   Md - distance matrix;
0015 %   lambda2 - algebraic connecity, i.e., lambda2(L)
0016 %   youtN - actual network size of resulting topology;
0017 %   L - the Laplacian matrix
0018 %
0019 % wzf, Dec 2008
0020 
0021 %   SynGrid
0022 %   Copyright (c) 2008, 2017-2018, Electric Power and Energy Systems (EPES) Research Lab
0023 %   by Zhifang Wang, Virginia Commonwealth University
0024 %
0025 %   This file is part of SynGrid.
0026 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0027 
0028 conn = 0; count =0;
0029 while(~conn) %-> this cycle for generation of connected toplogy
0030 
0031     % distance matrix Md computed by Md(i,j) = abs(i-j), the bus number distance
0032     Md=abs(((1:N)'*ones(1,N)-ones(N,1)*(1:N)));
0033     ds = nsw_get_nonzero(Md);
0034     switch N
0035         case 30
0036             d0 = 10;
0037         case 57
0038             d0 = 10;
0039         case 118
0040             d0 = 15;
0041         case 300
0042             d0 = 15; % 20, orignal setting, change it to 10 for bigger clustering coefficient
0043         otherwise
0044             d0 = 10;
0045     end
0046 
0047 
0048     % first step link selection -- among d0-neighbors
0049     Mda =or((Md>0 & Md<=d0),(Md>N/2 & Md>=(N-d0)));
0050     [Ma ]= nsw_clusterSW_sel_link2(Mda,degModel,mdeg); % choose links from local ngbrs
0051     degs = sum(or(Ma,Ma'),2);
0052     pos = find(degs==0);
0053 %     if(length(pos)>N*.1)
0054 %         continue;
0055 %     end
0056     % remove single isolated nodes
0057     [Ma, Md, Mda, Nn]= rmv_sgls_clusterSM(Ma,Md,Mda,N);
0058 
0059     rwd_yes = nsw_markov(alph,bta,Nn);
0060     rwdnodes_list = find(rwd_yes ==1);
0061     if(~isempty(rwdnodes_list))
0062         % decide rewiring clusters
0063         [clsts,clsts_size,n_clsts] = nsw_clusterSW_clstrs(rwd_yes,d0,Nn);
0064         rwd_clsts      = clsts;
0065 
0066         deg_all = sum(degs(rwdnodes_list));
0067         if(deg_all>0)
0068             p_sel = p_lrwd*mdeg*N*.5/deg_all;
0069             [Ma]= rewind(Ma,rwdnodes_list,rwd_clsts,p_sel);
0070         end
0071     end
0072     [conn,lambda2,L] = nsw_clusterSW_conn2(Ma);
0073     if(conn == 1)
0074        [A,indx,Ma] = Ma2net_clusterSM(or(Ma,Ma'),Md);
0075        youtN = Nn;
0076     end
0077 end
0078 
0079 %--------------------------------------------------------------------------
0080 %--------------------------------------------------------------------------
0081 function [Ma, Md, Mda,N]= rmv_sgls_clusterSM(Ma,Md,Mda,n)
0082 % remove the single isolated nodes from system
0083 
0084 degs = sum(or(Ma,Ma'),2);
0085 pos0 = find(degs==0); % nodes to remove
0086 for ii=1:length(pos0)
0087     ki=pos0(ii);
0088     ks_ngr = find(Mda(ki,:)>0);
0089     tmp =randi(length(ks_ngr));
0090     new_ngr = ks_ngr(tmp);
0091     Ma(ki,new_ngr)=1;
0092     Ma(new_ngr, ki)=1;
0093 end
0094 Ma = triu(Ma);
0095 %Md, Mda keep unchanged.
0096 N =n;
0097 % seq = 1:n; seq(pos0)=0;
0098 % pos = find(seq>0);   % nodes to keep
0099 % % shrink the vectors and matrix
0100 % tmp = Ma(:,pos); Ma = tmp(pos,:);
0101 % tmp = Md(:,pos); Md = tmp(pos,:);
0102 % tmp = Mda(:,pos); Mda = tmp(pos,:);
0103 % %x = x(pos); y = y(pos);
0104 % N = length(pos);
0105 
0106 %--------------------------------------------------------------------------
0107 %--------------------------------------------------------------------------
0108 function [Ma]= rewind(Ma,rwdnodes_list,rwd_clsts,p_sel)
0109 % rewind local links to across-network links
0110 
0111 Ma = or(Ma,Ma');
0112 n = length(Ma(:,1));
0113 for k = 1:length(rwdnodes_list)
0114     nodek = rwdnodes_list(k);
0115     ngsk  = find(Ma(nodek,:)>0);
0116     if(isempty(ngsk))
0117         continue;
0118     end
0119     degk = length(ngsk);
0120     pos  = find(rand(degk,1) < p_sel);
0121     if(~isempty(pos))
0122         old_ngsk = ngsk(pos);
0123         n_rwd = length(old_ngsk);
0124         new_ngsk = sel_newngs(rwd_clsts,nodek,n_rwd,n);
0125         
0126         if(length(new_ngsk)<n_rwd) % in case there are not enough rwd ngbrs
0127             n_rwd = length(new_ngsk);
0128             old_ngsk = old_ngsk(1:n_rwd);
0129         end
0130         
0131         
0132         Ma = update_Ma(Ma,nodek,old_ngsk,new_ngsk);
0133     end
0134 end
0135 
0136 
0137 %--------------------------------------------------------------------------
0138 function Ma_new = update_Ma(Ma,k,old_ngs,new_ngs)
0139 % update Ma, replace old neighbors with new neighbors
0140 
0141 Ma_new = Ma;
0142 
0143 % below section is added to avoid disconnectivity case
0144 old_ngK = sum(Ma(old_ngs,:),2); % node degree of old neighbors
0145 pos = find(old_ngK <= 1);       % if old ngbr only has this one link, not rewind it.
0146 if(~isempty(pos))
0147     old_ngs(pos)=0;
0148     new_pos = find(old_ngs>0);
0149     old_ngs = old_ngs(new_pos);   % shrink the rewiring list
0150     new_ngs = new_ngs(new_pos);   % shrink the rewiring list
0151 end
0152 Ma_new(k,old_ngs)= 0;
0153 Ma_new(k,new_ngs)=1;
0154 Ma_new(:,k)=Ma_new(k,:);
0155 
0156 %--------------------------------------------------------------------------
0157 function new_ngs = sel_newngs(rwd_clsts,nodek,n_rwd,n)
0158 % slect new neighbors for node-k from outside clsts
0159 
0160 sizec = size(rwd_clsts);
0161 nodes_avl =[];
0162 for ci = 1:sizec(1)
0163     nodes_ci = rwd_clsts{ci,1};
0164     if(isempty(find(nodes_ci == nodek)))
0165         nodes_avl = [nodes_avl nodes_ci];
0166     end
0167 end
0168 if(isempty(nodes_avl))
0169     nodes_aval = 1:n;
0170 end
0171 new_ngs = sel_nodes(nodes_avl,n_rwd);
0172 
0173 %--------------------------------------------------------------------------
0174 %--------------------------------------------------------------------------
0175 function nodes_sel = sel_nodes(nodes_avl,k)
0176 % select k nodes from nodes_aval
0177 m = length(nodes_avl);
0178 if(k<m)
0179     [tmp,ord] = sort(rand(1,m));
0180     nodes_sel = nodes_avl(ord(1:k));
0181 else
0182     nodes_sel = nodes_avl;
0183 end
0184 
0185 %--------------------------------------------------------------------------
0186 function yout = new2old_nodes(new_ord,nodes_new,clsts_new)
0187 % change new nodes number to old nodes number according to new_ord.
0188 %
0189 if(nargin == 2)
0190     m = length(nodes_new); nodes_old = zeros(1,m);
0191     for k = 1:m
0192         nodes_old(k) = find(new_ord == nodes_new(k));
0193     end
0194     yout = nodes_old;
0195 elseif(nargin == 3)
0196     sizenew= size(clsts_new);
0197     clsts_old = cell(sizenew);
0198     for ii = 1: sizenew(1)
0199         nodes_new = clsts_new{ii,1};
0200         nodes_old = new2old_nodes(new_ord,nodes_new);
0201         clsts_old{ii,1}= nodes_old;
0202     end
0203     yout = clsts_old;
0204 end
0205 
0206 %--------------------------------------------------------------------------
0207 function [A,indx,Ma] = Ma2net_clusterSM(Ma_full,Md)
0208 % Form the link admission matrix from the full connectivty matrix
0209 Ma = triu(Ma_full,1);
0210 N = length(Ma(:,1));
0211 [Lx,Ly] = find(Ma==1);
0212 num_lines = length(Lx);
0213 indx =[Lx Ly];
0214 % 4. form branch incidence matrix A
0215 A = sparse(zeros(num_lines,N));
0216 for ni = 1: num_lines
0217     A(ni, indx(ni,1)) = 1;
0218     A(ni, indx(ni,2)) = -1;
0219 end