fm_pmuloc.m

这是一个很适合研究和学习用的电力系统仿真软件

function [pmu_test, I_idx, pseudoi] = fm_pmuloc
% FM_PMULOC main function for PMU placement routines
%
% (PMU,I,PSEUDO_I) = FM_PMULOC
%      PMU placement set
%      I   measured currents
%      PSEUDO_I indirectly measured currents
%
% see also FM_PMUFIG
%
%Author:    Federico Milano
%Date:      11-Nov-2002
%Version:   1.0.0
%
%E-mail:    [e-mail]
%Web-site:  http://thunderbox.uwaterloo.ca/~fmilano
%
% Copyright (C) 2002-2006 Federico Milano

global Line Bus DAE Settings PMU
global Ltc File Varname Path Fig


if isempty(Line.con)
  fm_disp('No static network is loaded. LSSE is not intended for dynamic components.',1)
  return
end

tic

  if Fig.pmu
    hdl_pmv = findobj(Fig.pmu,'Tag','ListboxPMV');
    hdl_ang = findobj(Fig.pmu,'Tag','ListboxANG');
    hdl_V   = findobj(Fig.pmu,'Tag','ListboxV');
    hdl_pmc = findobj(Fig.pmu,'Tag','StaticTextPMC');
    hdl_pmu = findobj(Fig.pmu,'Tag','StaticTextPMU');
    hdl_nob = findobj(Fig.pmu,'Tag','StaticTextNOB');
    hdl_mv  = findobj(Fig.pmu,'Tag','StaticTextMV');
    hdl_mc  = findobj(Fig.pmu,'Tag','StaticTextMC');
  else
    hdl_pmv = 0;
    hdl_ang = 0;
    hdl_V   = 0;
    hdl_pmc = 0;
    hdl_pmu = 0;
    hdl_nob = 0;
    hdl_mv  = 0;
    hdl_mc  = 0;
  end

  type = PMU.method;

  fm_disp(' ',1)
  fm_disp('PMU Placement Routine',1)
  fm_disp(['Data File "',Path.data,File.data,'"'],1)

  P = DAE.glfp;
  Q = DAE.glfq;
  if Ltc.n > 0
    P(Bus.int(Ltc.con(:,1))) = P(Bus.int(Ltc.con(:,1))) - real(Ltc.dat(:,5));
    P(Bus.int(Ltc.con(:,2))) = P(Bus.int(Ltc.con(:,2))) - real(Ltc.dat(:,6));
    Q(Bus.int(Ltc.con(:,1))) = Q(Bus.int(Ltc.con(:,1))) + imag(Ltc.dat(:,5));
    Q(Bus.int(Ltc.con(:,2))) = Q(Bus.int(Ltc.con(:,2))) + imag(Ltc.dat(:,6));
  end
  P = round(abs(P)/Settings.lftol)*Settings.lftol;
  Q = round(abs(Q)/Settings.lftol)*Settings.lftol;
  zeroinj = P+Q;

  % vettore contenente i bus in cui si collocano i PMU
  pmu_con = [];
  % indice delle correnti misurate:
  % I_idx => [#corrente, #linea, from bus, to bus, sign]
  I_idx = [];

  % conteggio delle connessioni e ordinamento dei nodi
  nodi = [Line.con(:,1); Line.con(:,2)];
  n_link = zeros(Bus.n,2);
  for i = 1:Bus.n
    a = find(nodi == Bus.con(i,1));
    n_link(i,:) = [Bus.con(i,1), length(a)];
  end
  [y,i] = sort(n_link(:,2));
  n_link = n_link(i,:);
  linee = Line.con(:,[1 2]);
  pseudoi = 0;

  switch type

   case 1  % depth first method

    metodo = 'Depth First';
    fm_disp([metodo, ' Method '],1)

    while length(n_link(:,1)) > 0
      % locazione del PMU nel nodo non osservabile più interconnesso
      pmu_con = [pmu_con; n_link(length(n_link(:,1)),1)];
      i_idx = fm_iidx(pmu_con(length(pmu_con)),linee);
      I_idx = [I_idx; i_idx];
      % nodi osservabili dal primo PMU
      nodi_oss = [i_idx(:,4);pmu_con(length(pmu_con))];
      a = [];
      for i = 1:length(nodi_oss)
        a = [a; find(n_link(:,1) == nodi_oss(i))];
      end
      if ~isempty(a); n_link(a,:) = []; end
      if Fig.pmu
        set(hdl_pmu,'String',int2str(length(pmu_con)));
        set(hdl_nob,'String',int2str(length(n_link(:,1))));
        drawnow
      end
    end
    pmu_test = pmu_con;

   case 2  % depth first method with pseudo measurement of voltages and currents

    metodo = 'Graph Theoretic Procedure';
    fm_disp([metodo, ' Method '],1)

    P = DAE.glfp;
    Q = DAE.glfq;
    if Ltc.n > 0
      P(Bus.int(Ltc.con(:,1))) = P(Bus.int(Ltc.con(:,1))) - real(Ltc.dat(:,5));
      P(Bus.int(Ltc.con(:,2))) = P(Bus.int(Ltc.con(:,2))) - real(Ltc.dat(:,6));
      Q(Bus.int(Ltc.con(:,1))) = Q(Bus.int(Ltc.con(:,1))) + imag(Ltc.dat(:,5));
      Q(Bus.int(Ltc.con(:,2))) = Q(Bus.int(Ltc.con(:,2))) + imag(Ltc.dat(:,6));
    end
    P = abs(P);
    Q = abs(Q);
    nodi = [];
    while length(n_link(:,1)) > 0
      % locazione del PMU nel nodo non osservabile più interconnesso
      pmu_con = [pmu_con; n_link(length(n_link(:,1)),1)];
      i_idx = fm_iidx(pmu_con(length(pmu_con)),linee);

      if ~isempty(i_idx)
        I_idx = [I_idx; i_idx];
        % nodi osservabili dal primo PMU
        nodi_oss = [i_idx(:,4); pmu_con(length(pmu_con))];
        nodi = [nodi; nodi_oss];
        linee(i_idx(:,2),[1 2]) = zeros(length(nodi_oss)-1,2);
      end

      % determinazione delle pseudo-correnti determinate con la
      % legge di Kirchhoff per le correnti ed eliminazione dei nodi
      % di cui si può determinare la tensione con la legge di Ohm

      % determinazione delle pseudo-correnti nelle linee ai cui
      % estremi sono note le tensioni
      for ii = 1:length(nodi)
        I_idx_from = find(linee(:,1) == nodi(ii));
        I_idx_to = [];
        for jj = 1:length(nodi)
          ifrom = find(linee(I_idx_from,2) == nodi(jj));
          I_idx_to = [I_idx_to; I_idx_from(ifrom)];
        end
        if ~isempty(I_idx_to);
          n_current = length(I_idx_to);
          api = [[1:n_current]', I_idx_to];
          bpi = linee(I_idx_to,[1 2]);
          cpi = ones(length(I_idx_to),1);
          pi_idx = [api, bpi, cpi];
          linee(pi_idx(:,2),[1 2]) = zeros(length(pi_idx(:,1)),2);
          I_idx = [I_idx; pi_idx];
          pseudoi = pseudoi + length(pi_idx(:,1));
        end
      end

      count = 1;
      while count < length(nodi)

        if zeroinj(Bus.int(nodi(count))) == 0
          I_idx_from = find(linee(:,1) == nodi(count));
          I_idx_to = find(linee(:,2) == nodi(count));
          ncfrom = length(I_idx_from);
          ncto = length(I_idx_to);
          nc = ncfrom + ncto;
          if nc == 1
            if ncfrom == 1
              ki_idx = [1, I_idx_from, linee(I_idx_from,[1 2]), 1];
            else
              ki_idx = [1, I_idx_to, linee(I_idx_to,[2 1]), -1];
            end
            linee(ki_idx(2),[1 2]) = zeros(length(ki_idx(1)),2);
            I_idx = [I_idx; ki_idx];
            pseudoi = pseudoi + length(ki_idx(1));
            nodi_oss = [nodi_oss; ki_idx(4)];
            nodi = [ki_idx(4); nodi];

            % dopo l'aggiunta di un nuovo nodo misurato bisogna
            % ricontrollare tutti i nodi
            count = 1;

            % ricerca di pseudo-correnti che possono essere
            % misurate con il nodo aggiunto
            for ii = 1:length(nodi)
              I_idx_from = find(linee(:,1) == nodi(ii));
              I_idx_to = [];
              for jj = 1:length(nodi)
                ifrom = find(linee(I_idx_from,2) == nodi(jj));
                I_idx_to = [I_idx_to; I_idx_from(ifrom)];
              end
              if ~isempty(I_idx_to);
                n_current = length(I_idx_to);
                api = [[1:n_current]', I_idx_to];
                bpi = linee(I_idx_to,[1 2]);
                cpi = ones(length(I_idx_to),1);
                pi_idx = [api, bpi, cpi];
                linee(pi_idx(:,2),[1 2]) = zeros(length(pi_idx(:,1)),2);
                I_idx = [I_idx; pi_idx];
                pseudoi = pseudoi + length(pi_idx(:,1));
              end
            end

          else
            count = count + 1;
          end
        else
          count = count + 1;
        end
      end

      a = [];
      for i = 1:length(nodi_oss)
        a = [a; find(n_link(:,1) == nodi_oss(i))];
      end
      if ~isempty(a); n_link(a,:) = []; end
      if Fig.pmu
        set(hdl_pmu,'String',int2str(length(pmu_con)));
        set(hdl_nob,'String',int2str(length(n_link(:,1))));
        drawnow
      end
    end
    pmu_test = pmu_con;

   case 3  % annealing method

    metodo = 'Annealing Method';
    fm_disp([metodo, ' Method '],1)

    nodi = [];
    ok = 1;

    while length(n_link(:,1)) > 0
      % locazione del PMU nel nodo non osservabile più interconnesso
      pmu_con = [pmu_con; n_link(length(n_link(:,1)),1)];
      i_idx = fm_iidx(pmu_con(length(pmu_con)),linee);

      if ~isempty(i_idx)
        I_idx = [I_idx; i_idx];
        nodi_oss = [i_idx(:,4);pmu_con(length(pmu_con))];
        nodi = [nodi; nodi_oss];
        linee(i_idx(:,2),[1 2]) = zeros(length(nodi_oss)-1,2);
      end

      % determinazione delle pseudo-correnti  determinate con la
      % legge di Kirchhoff per le correnti ed eliminazione dei nodi
      % di cui si può determinare la tensione con la legge di Ohm

      % determinazione delle pseudo-correnti nelle linee ai cui
      % estremi sono note le tensioni
      for ii = 1:length(nodi)
        I_idx_from = find(linee(:,1) == nodi(ii));
        I_idx_to = [];
        for jj = 1:length(nodi)
          ifrom = find(linee(I_idx_from,2) == nodi(jj));
          I_idx_to = [I_idx_to; I_idx_from(ifrom)];
        end
        if ~isempty(I_idx_to);
          n_current = length(I_idx_to);
          api = [[1:n_current]', I_idx_to];
          bpi = linee(I_idx_to,[1 2]);
          cpi = ones(length(I_idx_to),1);
          pi_idx = [api, bpi, cpi];
          linee(pi_idx(:,2),[1 2]) = zeros(length(pi_idx(:,1)),2);
          I_idx = [I_idx; pi_idx];
          pseudoi = pseudoi + length(pi_idx(:,1));
        end
      end

      count = 1;
      while count < length(nodi)

        if zeroinj(Bus.int(nodi(count))) == 0
          I_idx_from = find(linee(:,1) == nodi(count));
          I_idx_to = find(linee(:,2) == nodi(count));
          ncfrom = length(I_idx_from);
          ncto = length(I_idx_to);
          nc = ncfrom + ncto;
          if nc == 1
            if ncfrom == 1
              ki_idx = [1, I_idx_from, linee(I_idx_from,[1 2]), 1];
            else
              ki_idx = [1, I_idx_to, linee(I_idx_to,[2 1]), -1];
            end
            linee(ki_idx(2),[1 2]) = zeros(length(ki_idx(1)),2);
            I_idx = [I_idx; ki_idx];
            pseudoi = pseudoi + length(ki_idx(1));
            nodi_oss = [nodi_oss; ki_idx(4)];
            nodi = [ki_idx(4); nodi];

            % dopo l'aggiunta di un nuovo nodo misurato bisogna
            % ricontrollare tutti i nodi
            count = 1;

            % ricerca di pseudo-correnti che possono essere
            % misurate con il nodo aggiunto
            for ii = 1:length(nodi)
              I_idx_from = find(linee(:,1) == nodi(ii));
              I_idx_to = [];
              for jj = 1:length(nodi)
                ifrom = find(linee(I_idx_from,2) == nodi(jj));
                I_idx_to = [I_idx_to; I_idx_from(ifrom)];
              end
              if ~isempty(I_idx_to);
                n_current = length(I_idx_to);
                api = [[1:n_current]', I_idx_to];
                bpi = linee(I_idx_to,[1 2]);
                cpi = ones(length(I_idx_to),1);
                pi_idx = [api, bpi, cpi];
                linee(pi_idx(:,2),[1 2]) = zeros(length(pi_idx(:,1)),2);
                I_idx = [I_idx; pi_idx];
                pseudoi = pseudoi + length(pi_idx(:,1));
              end
            end

          else
            count = count + 1;
          end
        else
          count = count + 1;
        end
      end

      a = [];
      for i = 1:length(nodi_oss)
        a = [a; find(n_link(:,1) == nodi_oss(i))];
      end
      if ~isempty(a); n_link(a,:) = []; end
      if Fig.pmu
        set(hdl_pmu,'String',int2str(length(pmu_con)));
        set(hdl_nob,'String',int2str(length(n_link(:,1))));
        drawnow
      end
    end

    pmu_test = pmu_con;
    pmu_test_old = pmu_con;
    I_idx_old = I_idx;
    pseudoi_old = pseudoi;
    nu = length(pmu_con);
    nl = 0;

    while (nu - nl) > 1 & nu > 1

      if nl == 0
        ntest = fix(0.85*nu);
      else
        ntest = fix((nu+nl)/2);
      end
      pmurand = randperm(nu);
      pmu_test = pmu_test_old(pmurand(1:ntest));

      % inizio della procedura di annealing
      linee = Line.con(:,[1 2]);
      I_idx = [];
      nodi = [];
      pseudoi = 0;

      for ii = 1:ntest
        i_idx = fm_iidx(pmu_test(ii),linee);