fm_connectivity.m

电力系统的psat

function [struct_connexe,connect_mat,bus_sub_net] = fm_connectivity(actual_time)

% function to detect islanded subnetwork after opening a breaker
% by Laurent Lenoir, June 2005.

%tic

global Line Bus Syn connect_mat struct_connexe nbr_sub_net
persistent old_time

% initialization
if actual_time == -999
  old_time = -999;
  return
end

% do not repeat computations if the simulation is stucking
if old_time ~= actual_time
  old_time = actual_time;
else
  return
end

struct_connexe = cell(1,Bus.n);
bus_sub_net = {};
bus_disconnected = [];
bus_discon_mat = [];

% connectivity matrix
connect_mat = sparse(Line.from,Line.to,1,Bus.n,Bus.n);
connect_mat = connect_mat + sparse(Line.to,Line.from,1,Bus.n,Bus.n);

% check for open lines
line_discon = find(Line.con(:,9) > 1e5);

if ~isempty(line_discon)

  % update connectivity matrix
  for ik = 1:length(line_discon)
    hhh = Line.from(line_discon(ik));
    kkk = Line.to(line_discon(ik));
    connect_mat(hhh,kkk) = connect_mat(hhh,kkk) - 1;
    connect_mat(kkk,hhh) = connect_mat(kkk,hhh) - 1;
  end

  % search link between sub_network
  for ik = 1:length(line_discon)
    SommetsAccessibles(Bus.n,Line.to(line_discon(ik)));
    SommetsAccessibles(Bus.n,Line.from(line_discon(ik)));
    bus_disconnected = [bus_disconnected; ...
                        Line.from(line_discon(ik)); ...
                        Line.to(line_discon(ik))];
  end
  
  % create a matrix ij of bus deconnected with a link
  bus_discon_mat = zeros(length(bus_disconnected),length(bus_disconnected));
  for jk = 1:length(bus_disconnected)
    gh = zeros(length(bus_disconnected),1);
    for ik = 1:length(bus_disconnected)
      hg = find(struct_connexe{bus_disconnected(jk)}==bus_disconnected(ik));
      if ~isempty(hg)
        gh(ik) = bus_disconnected(ik);
      end
    end
    bus_discon_mat(:,jk) = gh;
  end
  
  % elimination sous reseau identique
  jk = 1;
  while jk <= length(bus_discon_mat(1,:))
    [ih,jh,vh] = find(bus_discon_mat(:,jk));
    if length(jh)>1
      [ih1,jh1]=find(bus_discon_mat==vh(1));
      for kl=2:length(jh1)
        %bus_discon_mat = delete_column(bus_discon_mat,jh1(kl));
        bus_discon_mat(:,jh1(kl)) = [];
      end
    end
    jk = jk + 1;
  end
  
  nbr_sub_net=length(bus_discon_mat(1,:));
  
  if (nbr_sub_net > 1)
    disp(['There are ' num2str(nbr_sub_net) ' islanded network'])
  end
  
  for ik=1:nbr_sub_net
    [ih,jh,vh]=find(bus_discon_mat(:,ik));
    bus_sub_net{ik}=struct_connexe{vh(1)};
    nbr_gen_sub_network=proc_nbr_gen(bus_sub_net{ik});
    if nbr_sub_net > 1
      disp(['There is a sub_network with ', ...
            num2str(length(bus_sub_net{ik})), ...
            ' buses with ', num2str(nbr_gen_sub_network), ...
            ' generator(s)'])
    end
  end
  
%toc
end

%--------------------------------------------------------------------------
function Visite(IndexSommet)

global ComptrVisite NumeroDOrdre NbSommets connect_mat struct_connexe lk

ComptrVisite=ComptrVisite+1;
NumeroDOrdre(IndexSommet)=ComptrVisite;
struct_connexe{lk}(ComptrVisite,1)=IndexSommet;

for AutreSommet=1:NbSommets
  if ((connect_mat(IndexSommet,AutreSommet)>0) & ...
      (NumeroDOrdre(AutreSommet)==0))
    Visite(AutreSommet);
  end
end

%--------------------------------------------------------------------------
function SommetsAccessibles(NbSommet,Sommet)

global ComptrVisite NbSommets NumeroDOrdre connect_mat struct_connexe lk

NumeroDOrdre=zeros(length(connect_mat(:,1)),1);
NbSommets=NbSommet;
lk=Sommet;
ComptrVisite=0;
Visite(Sommet);

%--------------------------------------------------------------------------
function nbr_gen = proc_nbr_gen(bus_sub_network)

global Syn

nbr_gen=0;
for ik=1:length(bus_sub_network)
  hg=find(Syn.bus==bus_sub_network(ik));
  if ~isempty(hg)
    nbr_gen=nbr_gen+1;
  end
end